#include "Mesh.h"
#include "IO/IO.h"

//#include "cgnslib.h"

#include <math.h>
#include <set>


Mesh::Mesh(Config* config_) :
    config(config_)
{
    
    //readGmshMesh();

    //一维、二维、三维网格，均由这个一个函数读取
    readCgnsMesh();

    if (config->dimension == 3)
    {
        setInternalFaceGroup();

        setMaterial();

        //creating mfem mesh, mfem periodic mesh
        init();

        calConnectivity();

        //log_info << "end creating serial mesh" << std::endl;
    }
    else if (config->dimension == 2)
    {
        setInternalLineGroup();

        setMaterial2d();

        //creating mfem mesh, mfem periodic mesh
        if (config->debugLevel > 1) log_info << "init2d begin";
        init2d();
        if (config->debugLevel > 1) log_info << "init2d end";

        if (config->debugLevel > 1) log_info << "calConnectivity2d begin";
        calConnectivity2d();
        if (config->debugLevel > 1) log_info << "calConnectivity2d end";

        //log_info << "end creating serial mesh" << std::endl;
    }
    else
    {
        log_error << "not surppot dimension: " << config->dimension;
        exit(0);
    }


}




void Mesh::readCgnsMesh()
{
    //if(config->debugLevel > 0) log_info << "---------- begin: read cgns mesh";

    ////cgns的单元是一起存储的，包括line face volume，需要将line face volume分开存储，重新编号，因为EasyPIC中faceGroup中是face的编号
    //int index_file,index_base,index_zone,nbocos,ib;

    //cgsize_t isize[3][1],ielem[20*16*8][8];

    //int normalindex[3],ndataset;
    //int i,normallist;

    //char boconame[33];

    //cgsize_t irmin,irmax,istart,iend;

    //int nsections,index_sect,nbndry,iparent_flag;
    //cgsize_t iparentdata;

    //char zonename[33],sectionname[33];

    //CGNS_ENUMT(ElementType_t) itype;

    //CGNS_ENUMT(BCType_t) ibocotype;
    //CGNS_ENUMT(PointSetType_t) iptset;
    //CGNS_ENUMT(DataType_t) normaldatatype;
    //CGNS_ENUMT(GridLocation_t) igr;
    //cgsize_t npts,normallistflag;



    ////open CGNS file for read-only
    //if (cg_open(config->mesh_param.meshName.c_str(),CG_MODE_READ,&index_file)) cg_error_exit();

    ////we know there is only one base (real working code would check!)
    //index_base=1;
    //// we know there is only one zone (real working code would check!)
    //index_zone=1;

    //// get zone size (and name - although not needed here)
    //cg_zone_read(index_file,index_base,index_zone,zonename,isize[0]);

    //// lower range index
    //irmin=1;
    //// upper range index of vertices
    //irmax=isize[0][0];
    //std::vector<double> x,y,z;
    //x.resize(irmax);
    //y.resize(irmax);
    //z.resize(irmax);

    //// read grid coordinates
    //cg_coord_read(index_file,index_base,index_zone,"CoordinateX",
    //              CGNS_ENUMV(RealDouble),&irmin,&irmax,&x[0]);
    //cg_coord_read(index_file,index_base,index_zone,"CoordinateY",
    //              CGNS_ENUMV(RealDouble),&irmin,&irmax,&y[0]);
    //cg_coord_read(index_file,index_base,index_zone,"CoordinateZ",
    //              CGNS_ENUMV(RealDouble),&irmin,&irmax,&z[0]);

    //pointList.resize(x.size());
    //for(int iPoint = 0; iPoint < pointList.size(); iPoint++)
    //{
    //    Point& point = pointList[iPoint];
    //    point.position(0) = x[iPoint] * config->mesh_param.meshScale;
    //    point.position(1) = y[iPoint] * config->mesh_param.meshScale;
    //    point.position(2) = z[iPoint] * config->mesh_param.meshScale;
    //}



    ////find out how many sections
    //cg_nsections(index_file,index_base,index_zone,&nsections);
    ////log_info << "number of sections = " << nsections;

    //
    ////cgns网格的faceGroup等存储的所有单元的序号，而EasyPSim是存储某类单元的序号，比如LineGroup存储的是LineList的序号
    ////faceGroup存储的是FaceList中序号，因此需要新建临时map存储对应关系
    //std::map<int, int> lineIndexGlobalToLocalMap;
    //std::map<int, int> faceIndexGlobalToLocalMap;
    //std::map<int, int> volumeIndexGlobalToLocalMap;


    //// read element connectivity
    //int iElement = 0;
    //for (index_sect=1; index_sect <= nsections; index_sect++)
    //{
    //    cg_section_read(index_file,index_base,index_zone,index_sect,sectionname,
    //                    &itype,&istart,&iend,&nbndry,&iparent_flag);
    //    //log_info << "Reading section data...";
    //    //log_info << "section name = ",sectionname;
    //    //log_info << "section type = ", ElementTypeName[itype];
    //    //log_info << "istart,iend = " << (int)istart << " " << (int)iend;

    //    if (itype == CGNS_ENUMV(HEXA_8))
    //    {
    //        log_error << "not support HEXA_8 element";
    //        exit(0);
    //    }
    //    else if (itype == CGNS_ENUMV(TETRA_4))
    //    {
    //        int volumeNum = (int)iend - (int)istart + 1;
    //        cgsize_t* pointIndex = new cgsize_t[volumeNum * 4];

    //        cg_elements_read(index_file,index_base,index_zone,index_sect,pointIndex,&iparentdata);

    //        for(int iVolume = 0; iVolume < volumeNum; iVolume++)
    //        {
    //            std::vector<int> elementListTemp;
    //            elementListTemp.resize(4);
    //            elementListTemp[0] = pointIndex[iVolume * 4 + 0] - 1;
    //            elementListTemp[1] = pointIndex[iVolume * 4 + 1] - 1;
    //            elementListTemp[2] = pointIndex[iVolume * 4 + 2] - 1;
    //            elementListTemp[3] = pointIndex[iVolume * 4 + 3] - 1;
    //            Element* element = new Element(geometry::TETRAHEDRON, elementListTemp);
    //            element->elementIndex = iElement;
    //            elementList.push_back(element);
    //            iElement++;


    //            //cgns index start from 1
    //            Volume volumeTemp;
    //            volumeTemp.pointList[0] = pointIndex[iVolume*4+0] - 1;
    //            volumeTemp.pointList[1] = pointIndex[iVolume*4+1] - 1;
    //            volumeTemp.pointList[2] = pointIndex[iVolume*4+2] - 1;
    //            volumeTemp.pointList[3] = pointIndex[iVolume*4+3] - 1;
    //            volumeList.push_back(volumeTemp);

    //            volumeIndexGlobalToLocalMap[iElement - 1] = volumeList.size() - 1;

    //        }

    //        delete[] pointIndex;


    //    }
    //    else if (itype == CGNS_ENUMV(TRI_3))
    //    {
    //        //log_info << "this is face: TRI_3";

    //        int faceNum = (int)iend - (int)istart + 1;
    //        cgsize_t* pointIndex = new cgsize_t[faceNum * 3];

    //        cg_elements_read(index_file,index_base,index_zone,index_sect,pointIndex,&iparentdata);

    //        for(int iFace = 0; iFace < faceNum; iFace++)
    //        {
    //            std::vector<int> elementListTemp;
    //            elementListTemp.resize(3);
    //            elementListTemp[0] = pointIndex[iFace * 3 + 0] - 1;
    //            elementListTemp[1] = pointIndex[iFace * 3 + 1] - 1;
    //            elementListTemp[2] = pointIndex[iFace * 3 + 2] - 1;
    //            Element* element = new Element(geometry::TRIANGLE, elementListTemp);
    //            element->elementIndex = iElement;
    //            elementList.push_back(element);
    //            iElement++;

    //            //cgns index start from 1
    //            Face faceTemp;
    //            faceTemp.pointList[0] = pointIndex[iFace * 3 + 0] - 1;
    //            faceTemp.pointList[1] = pointIndex[iFace * 3 + 1] - 1;
    //            faceTemp.pointList[2] = pointIndex[iFace * 3 + 2] - 1;
    //            faceList.push_back(faceTemp);
    //            
    //            faceIndexGlobalToLocalMap[iElement - 1] = faceList.size() - 1;

    //        }

    //        delete[] pointIndex;
    //    }
    //    else if (itype == CGNS_ENUMV(BAR_2))
    //    {
    //        int lineNum = (int)iend - (int)istart + 1;
    //        cgsize_t* pointIndex = new cgsize_t[lineNum * 2];

    //        cg_elements_read(index_file, index_base, index_zone, index_sect, pointIndex, &iparentdata);

    //        for (int iLine = 0; iLine <  lineNum; iLine++)
    //        {
    //            std::vector<int> elementListTemp;
    //            elementListTemp.resize(3);
    //            elementListTemp[0] = pointIndex[iLine * 2 + 0] - 1;
    //            elementListTemp[1] = pointIndex[iLine * 2 + 1] - 1;
    //            Element* element = new Element(geometry::SEGMENT, elementListTemp);
    //            element->elementIndex = iElement;
    //            elementList.push_back(element);
    //            iElement++;

    //            //cgns index start from 1
    //            Line lineTemp;
    //            lineTemp.pointList[0] = pointIndex[iLine * 2 + 0] - 1;
    //            lineTemp.pointList[1] = pointIndex[iLine * 2 + 1] - 1;
    //            lineList.push_back(lineTemp);

    //            lineIndexGlobalToLocalMap[iElement - 1] = lineList.size() - 1;

    //        }

    //        delete[] pointIndex;
    //    }
    //    else
    //    {
    //        log_error << "not support such type element: " << itype;
    //        exit(0);
    //    }
    //}


    //// find out number of BCs that exist under this zone
    //cg_nbocos(index_file,index_base,index_zone,&nbocos);

    ////log_info << "reading mesh boundary, number of boundary = " << nbocos;
    ////do loop over the total number of BCs
    //for (ib=1; ib <= nbocos; ib++)
    //{
    //    //find out what BC grid location is (expecting FaceCenter)
    //    cg_goto(index_file,index_base,"Zone_t",1,"ZoneBC_t",1,"BC_t",ib,"end");
    //    cg_gridlocation_read(&igr);
    //    if (igr == CGNS_ENUMV(FaceCenter))
    //    {
    //        if (config->dimension == 3)
    //        {
    //            //FaceCenter：表示这是三维网格的边界faceGroup
    //            if (config->debugLevel > 1) log_info << "this is faceGroup, igr = " << igr;

    //            std::string faceGroupName;

    //            //get BC info
    //            cg_boco_info(index_file, index_base, index_zone, ib, boconame, &ibocotype,
    //                &iptset, &npts, normalindex, &normallistflag, &normaldatatype, &ndataset);
    //            if (iptset != CGNS_ENUMV(PointList))
    //            {
    //                log_info << "Error.  For this program, BCs must be set up as PointList type: " << PointSetTypeName[iptset];
    //            }
    //            //log_info << "BC number: " << ib;
    //            //log_info << "   name = " << boconame;
    //            //log_info << "   type = " << BCTypeName[ibocotype];
    //            //log_info << "   no of elements = " << (int)npts;

    //            faceGroupName = boconame;
    //            cgsize_t* faceIndex = new cgsize_t[(int)npts];

    //            /* read point list in here (nothing done with them in this program) */
    //            cg_boco_read(index_file, index_base, index_zone, ib, faceIndex, &normallist);
    //            //log_info << "      these elements read here, but only some printed out: ";

    //            FaceGroup faceGroup;
    //            faceGroup.faceList.resize((int)npts);
    //            faceGroup.name = faceGroupName;

    //            for (int i = 0; i < (int)npts; i++)
    //            {
    //                //cgns index start from 1
    //                faceGroup.faceList[i] = faceIndexGlobalToLocalMap.find(faceIndex[i] - 1)->second;
    //                faceGroup.faceElementList.push_back(elementList[faceIndex[i] - 1]);
    //            }

    //            log_info << "faceGroup: " << faceGroupName;
    //            faceGroupMap[faceGroupName] = faceGroup;

    //            delete[] faceIndex;
    //        }
    //        else if (config->dimension == 2)
    //        {
    //            //EdgeCenter: 表示这是二维网格的边界lineGroup
    //            if (config->debugLevel > 1) log_info << "this is lineGroup, igr = " << igr;

    //            std::string lineGroupName;

    //            //get BC info
    //            cg_boco_info(index_file, index_base, index_zone, ib, boconame, &ibocotype,
    //                &iptset, &npts, normalindex, &normallistflag, &normaldatatype, &ndataset);
    //            if (iptset != CGNS_ENUMV(PointList))
    //            {
    //                log_info << "Error.  For this program, BCs must be set up as PointList type: " << PointSetTypeName[iptset];
    //            }
    //            //log_info << "BC number: " << ib;
    //            //log_info << "   name = " << boconame;
    //            //log_info << "   type = " << BCTypeName[ibocotype];
    //            //log_info << "   no of elements = " << (int)npts;

    //            lineGroupName = boconame;
    //            cgsize_t* lineIndex = new cgsize_t[(int)npts];

    //            /* read point list in here (nothing done with them in this program) */
    //            cg_boco_read(index_file, index_base, index_zone, ib, lineIndex, &normallist);
    //            //log_info << "      these elements read here, but only some printed out: ";

    //            LineGroup lineGroup;
    //            lineGroup.lineList.resize((int)npts);
    //            lineGroup.name = lineGroupName;
    //            for (int i = 0; i < (int)npts; i++)
    //            {
    //                //cgns index start from 1
    //                lineGroup.lineList[i] = lineIndexGlobalToLocalMap.find(lineIndex[i] - 1)->second;;
    //                lineGroup.lineElementList.push_back(elementList[lineIndex[i] - 1]);
    //            }

    //            //get volumeGroup material
    //            //std::string materialName = config->mesh_param.volumeGroupToMaterialMap[lineGroup.name];
    //            //lineGroup.material = config->materialMap[materialName];
    //            log_info << "lineGroup: " << lineGroupName;
    //            lineGroupMap[lineGroupName] = lineGroup;

    //            delete[] lineIndex;
    //        }
    //        else
    //        {
    //            log_error << "dimension error, " << config->dimension;

    //        }
    //    }
    //    else if(igr == CGNS_ENUMV(CellCenter))
    //    {
    //        if (config->dimension == 3)
    //        {
    //            //CellCenter: 表示这是三维网格的子网格volumeGroup，比如介质球内部网格
    //            if (config->debugLevel > 1) log_info << "this is volumeGroup, igr = " << igr;


    //            std::string volumeGroupName;

    //            //get BC info
    //            cg_boco_info(index_file, index_base, index_zone, ib, boconame, &ibocotype,
    //                &iptset, &npts, normalindex, &normallistflag, &normaldatatype, &ndataset);
    //            if (iptset != CGNS_ENUMV(PointList))
    //            {
    //                log_info << "Error.  For this program, BCs must be set up as PointList type: " << PointSetTypeName[iptset];
    //            }
    //            //log_info << "BC number: " << ib;
    //            //log_info << "   name = " << boconame;
    //            //log_info << "   type = " << BCTypeName[ibocotype];
    //            //log_info << "   no of elements = " << (int)npts;


    //            volumeGroupName = boconame;
    //            cgsize_t* volumeIndex = new cgsize_t[(int)npts];

    //            /* read point list in here (nothing done with them in this program) */
    //            cg_boco_read(index_file, index_base, index_zone, ib, volumeIndex, &normallist);

    //            VolumeGroup volumeGroup;
    //            volumeGroup.volumeList.resize((int)npts);
    //            volumeGroup.name = volumeGroupName;
    //            for (int i = 0; i < (int)npts; i++)
    //            {
    //                //cgns index start from 1
    //                volumeGroup.volumeList[i] = volumeIndexGlobalToLocalMap.find(volumeIndex[i] - 1)->second;;
    //                volumeGroup.volumeElementList.push_back(elementList[volumeIndex[i] - 1]);
    //            }

    //            log_info << "volumeGroup: " << volumeGroupName;
    //            volumeGroupMap[volumeGroupName] = volumeGroup;

    //            delete[] volumeIndex;
    //        }
    //        else if (config->dimension == 2)
    //        {
    //            //FaceCenter：表示这是三维网格的边界faceGroup
    //            if (config->debugLevel > 1) log_info << "this is faceGroup, igr = " << igr;

    //            std::string faceGroupName;

    //            //get BC info
    //            cg_boco_info(index_file, index_base, index_zone, ib, boconame, &ibocotype,
    //                &iptset, &npts, normalindex, &normallistflag, &normaldatatype, &ndataset);
    //            if (iptset != CGNS_ENUMV(PointList))
    //            {
    //                log_info << "Error.  For this program, BCs must be set up as PointList type: " << PointSetTypeName[iptset];
    //            }
    //            //log_info << "BC number: " << ib;
    //            //log_info << "   name = " << boconame;
    //            //log_info << "   type = " << BCTypeName[ibocotype];
    //            //log_info << "   no of elements = " << (int)npts;

    //            faceGroupName = boconame;
    //            cgsize_t* faceIndex = new cgsize_t[(int)npts];

    //            /* read point list in here (nothing done with them in this program) */
    //            cg_boco_read(index_file, index_base, index_zone, ib, faceIndex, &normallist);
    //            //log_info << "      these elements read here, but only some printed out: ";

    //            FaceGroup faceGroup;
    //            faceGroup.faceList.resize((int)npts);
    //            faceGroup.name = faceGroupName;

    //            for (int i = 0; i < (int)npts; i++)
    //            {
    //                //cgns index start from 1
    //                faceGroup.faceList[i] = faceIndexGlobalToLocalMap.find(faceIndex[i] - 1)->second;
    //                faceGroup.faceElementList.push_back(elementList[faceIndex[i] - 1]);
    //            }

    //            log_info << "faceGroup: " << faceGroupName;
    //            faceGroupMap[faceGroupName] = faceGroup;

    //            delete[] faceIndex;
    //        }
    //        else
    //        {
    //            log_error << "dimension error, " << config->dimension;
    //        }

    //    }
    //    else if (igr == CGNS_ENUMV(EdgeCenter))
    //    {
    //        //EdgeCenter: 表示这是二维网格的边界lineGroup
    //        if(config->debugLevel > 1) log_info << "this is lineGroup, igr = " << igr;

    //        std::string lineGroupName;

    //        //get BC info
    //        cg_boco_info(index_file, index_base, index_zone, ib, boconame, &ibocotype,
    //            &iptset, &npts, normalindex, &normallistflag, &normaldatatype, &ndataset);
    //        if (iptset != CGNS_ENUMV(PointList))
    //        {
    //            log_info << "Error.  For this program, BCs must be set up as PointList type: " << PointSetTypeName[iptset];
    //        }
    //        //log_info << "BC number: " << ib;
    //        //log_info << "   name = " << boconame;
    //        //log_info << "   type = " << BCTypeName[ibocotype];
    //        //log_info << "   no of elements = " << (int)npts;

    //        lineGroupName = boconame;
    //        cgsize_t* lineIndex = new cgsize_t[(int)npts];

    //        /* read point list in here (nothing done with them in this program) */
    //        cg_boco_read(index_file, index_base, index_zone, ib, lineIndex, &normallist);
    //        //log_info << "      these elements read here, but only some printed out: ";

    //        LineGroup lineGroup;
    //        lineGroup.lineList.resize((int)npts);
    //        lineGroup.name = lineGroupName;
    //        for (int i = 0; i < (int)npts; i++)
    //        {
    //            //cgns index start from 1
    //            lineGroup.lineList[i] = lineIndexGlobalToLocalMap.find(lineIndex[i] - 1)->second;;
    //            lineGroup.lineElementList.push_back(elementList[lineIndex[i] - 1]);
    //        }

    //        //get volumeGroup material
    //        //std::string materialName = config->mesh_param.volumeGroupToMaterialMap[lineGroup.name];
    //        //lineGroup.material = config->materialMap[materialName];

    //        log_info << "lineGroup: " << lineGroupName;
    //        lineGroupMap[lineGroupName] = lineGroup;

    //        delete[] lineIndex;
    //    }


    //}


    if (config->debugLevel > 0) log_info << "---------- end: read cgns mesh";

}


void Mesh::init()
{
    
    //============================== begin: creating serial mfem mesh ======================
    if(config->debugLevel > 0) log_info << "begin: creating mfem serial mesh";

    //mpi init at ParallelMPI.cpp
    //mfem::MPI_Session mfem::mpi;

    const char* device_config = "cpu";
    mfem::Device device(device_config);

    double* vertices;
    int num_vertices = pointList.size();
    vertices = new double[num_vertices * 3];


    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        auto& point = pointList[iPoint];
        vertices[3 * iPoint + 0] = point.position(0);
        vertices[3 * iPoint + 1] = point.position(1);
        vertices[3 * iPoint + 2] = point.position(2);
    }

    std::vector<int> element_indices;
    std::vector<int> element_attributes;
    mfem::Geometry::Type element_type = mfem::Geometry::TETRAHEDRON;
    int num_elements = volumeList.size();

    element_indices.resize(volumeList.size() * 4);
    element_attributes.resize(volumeList.size());
    for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        element_attributes[iVolume] = 0;
        element_indices[iVolume * 4 + 0] = volumeList[iVolume].pointList[0];
        element_indices[iVolume * 4 + 1] = volumeList[iVolume].pointList[1];
        element_indices[iVolume * 4 + 2] = volumeList[iVolume].pointList[2];
        element_indices[iVolume * 4 + 3] = volumeList[iVolume].pointList[3];

        element_attributes[iVolume] = 1;
    }
    //** attribute for element
    //cout << "volumeGroupMap.size-----------------------------------" << volumeGroupMap.size() << endl;
    int index = 10;
    for (auto iter = volumeGroupMap.begin(); iter != volumeGroupMap.end(); iter++)
    {
        auto volumeGroupName = iter->first;
        auto& volumeGroup = iter->second;
        volumeGroup.attributeIndex = index;
        index = index + 1;
        for (int i = 0; i < volumeGroup.volumeList.size(); i++)
        {
            //cout << "volumeGroup.volumeList[i]======="<<i<<"======" << volumeGroup.volumeList[i] << endl;
            element_attributes[volumeGroup.volumeList[i]] = volumeGroup.attributeIndex;
        }
    }



    std::vector<int> boundary_indices;
    std::vector<int> boundary_attributes;
    int num_boundary_elements = 0;
    mfem::Geometry::Type boundary_type = mfem::Geometry::TRIANGLE;

    int attributeIndex = 1;
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        faceGroup.attributeIndex = attributeIndex;
        if (config->debugLevel > 1) log_info << faceGroupName << " attribute : " << faceGroup.attributeIndex;

        attributeIndex++;

        //mfem not support internal face
        if (faceGroup.isInternal && config->electromagneticFieldType == "electrostatic")
        {
            continue;
        }

        auto& faceElementList = faceGroup.faceElementList;
        for (int iFaceElement = 0; iFaceElement < faceElementList.size(); iFaceElement++)
        {
            Element* faceElement = faceElementList[iFaceElement];
            boundary_indices.push_back(faceElement->pointList[0]);
            boundary_indices.push_back(faceElement->pointList[1]);
            boundary_indices.push_back(faceElement->pointList[2]);

            boundary_attributes.push_back(faceGroup.attributeIndex);
        }
        
    }
    num_boundary_elements = boundary_attributes.size();


    int dimension = config->dimension;
    int order = config->order;


    
    //============================== end: creating serial mfem mesh ======================



    // calculate faceList from mfem::Mesh
    int numFace = mfemMesh->GetNFaces();
    //log_info << "serial mesh (maybe periodic) vertex number: " << mfemMesh->GetNV();
    //log_info << "serial mesh (maybe periodic) face number: " << numFace;
    //log_info << "serial mesh (maybe periodic) volume number: " << mfemMesh->GetNE();
    faceList.resize(numFace);
    for (int iFace = 0; iFace < numFace; iFace++)
    {
        const mfem::Element* element = mfemMesh->GetFace(iFace);
        const int* vertexIndexs;
        vertexIndexs = element->GetVertices();
        faceList[iFace].pointList[0] = vertexIndexs[0];
        faceList[iFace].pointList[1] = vertexIndexs[1];
        faceList[iFace].pointList[2] = vertexIndexs[2];
    }

    numberPoint = pointList.size();
    numberLine = lineList.size();
    numberFace = faceList.size();
    numberVolume = volumeList.size();


    //cal maxGeometryLength,maxGeometryPosition
    maxGeometryLength = 0.0;
    minGeometryPoint.position(0) = 1.0e99;
    minGeometryPoint.position(1) = 1.0e99;
    minGeometryPoint.position(2) = 1.0e99;
    maxGeometryPoint.position(0) = -1.0e99;
    maxGeometryPoint.position(1) = -1.0e99;
    maxGeometryPoint.position(2) = -1.0e99;

    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        Point& point = pointList[iPoint];
        if (point.position(0) < minGeometryPoint.position(0))
        {
            minGeometryPoint.position(0) = point.position(0);
        }
        if (point.position(1) < minGeometryPoint.position(1))
        {
            minGeometryPoint.position(1) = point.position(1);
        }
        if (point.position(2) < minGeometryPoint.position(2))
        {
            minGeometryPoint.position(2) = point.position(2);
        }

        if (point.position(0) > maxGeometryPoint.position(0))
        {
            maxGeometryPoint.position(0) = point.position(0);
        }
        if (point.position(1) > maxGeometryPoint.position(1))
        {
            maxGeometryPoint.position(1) = point.position(1);
        }
        if (point.position(2) > maxGeometryPoint.position(2))
        {
            maxGeometryPoint.position(2) = point.position(2);
        }
    }


    maxGeometryLength = calLineLength(minGeometryPoint, maxGeometryPoint);

    //originOfCoordinate.position(0) = 0.5 * (minGeometryPoint.position(0) + maxGeometryPoint.position(0));
    //originOfCoordinate.position(1) = 0.5 * (minGeometryPoint.position(1) + maxGeometryPoint.position(1));
    //originOfCoordinate.position(2) = 0.5 * (minGeometryPoint.position(2) + maxGeometryPoint.position(2));
    radius = 0.5 * (maxGeometryPoint.position(0) - minGeometryPoint.position(0));

    //calculate minLineLength maxLineLength
    minLineLength = maxGeometryLength;
    maxLineLength = 0.0;
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        Face& face = faceList[iFace];
        for (int i = 0; i < 3; i++)
        {
            Point& point0 = pointList[face.pointList[(i + 1) % 3]];
            Point& point1 = pointList[face.pointList[(i + 2) % 3]];
            double length = calLineLength(point0, point1);
            if (length < minLineLength)
            {
                minLineLength = length;
            }
            if (length > maxLineLength)
            {
                maxLineLength = length;
            }
        }
    }
    //log_info << "minLineLength: " << minLineLength << " " << faceList.size();

    lengthSmallValue = smallValue * minLineLength;
    lengthVerySmallValue = verySmallValue * minLineLength;
    lengthLargeValue = largeValue * maxLineLength;
    lengthVeryLargeValue = veryLargeValue * maxLineLength;




    if (config->dimension == 3)
    {
        numberCell = numberVolume;
    }
    else
    {
        log_error << "dimension no right: " << config->dimension;
        exit(0);
    }
    if (config->debugLevel > 1) log_info << "end: init";

}


void Mesh::init2d()
{
    
    //cal maxGeometryLength,maxGeometryPosition
    maxGeometryLength = 0.0;
    minGeometryPoint.position(0) = 1.0e99;
    minGeometryPoint.position(1) = 1.0e99;
    minGeometryPoint.position(2) = 1.0e99;
    maxGeometryPoint.position(0) = -1.0e99;
    maxGeometryPoint.position(1) = -1.0e99;
    maxGeometryPoint.position(2) = -1.0e99;

    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        Point& point = pointList[iPoint];
        if (point.position(0) < minGeometryPoint.position(0))
        {
            minGeometryPoint.position(0) = point.position(0);
        }
        if (point.position(1) < minGeometryPoint.position(1))
        {
            minGeometryPoint.position(1) = point.position(1);
        }
        if (point.position(2) < minGeometryPoint.position(2))
        {
            minGeometryPoint.position(2) = point.position(2);
        }

        if (point.position(0) > maxGeometryPoint.position(0))
        {
            maxGeometryPoint.position(0) = point.position(0);
        }
        if (point.position(1) > maxGeometryPoint.position(1))
        {
            maxGeometryPoint.position(1) = point.position(1);
        }
        if (point.position(2) > maxGeometryPoint.position(2))
        {
            maxGeometryPoint.position(2) = point.position(2);
        }
    }
    maxGeometryLength = calLineLength(minGeometryPoint, maxGeometryPoint);

    //originOfCoordinate.position(0) = 0.5 * (minGeometryPoint.position(0) + maxGeometryPoint.position(0));
    //originOfCoordinate.position(1) = 0.5 * (minGeometryPoint.position(1) + maxGeometryPoint.position(1));
    //originOfCoordinate.position(2) = 0.5 * (minGeometryPoint.position(2) + maxGeometryPoint.position(2));
    radius = 0.5 * (maxGeometryPoint.position(0) - minGeometryPoint.position(0));

    if (config->debugLevel > 1) log_info << "maxGeometryLength = " << maxGeometryLength;


    //calculate minLineLength maxLineLength
    minLineLength = maxGeometryLength;
    maxLineLength = 0.0;
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        Face& face = faceList[iFace];
        for (int i = 0; i < 3; i++)
        {
            Point& point0 = pointList[face.pointList[(i + 1) % 3]];
            Point& point1 = pointList[face.pointList[(i + 2) % 3]];
            double length = calLineLength(point0, point1);
            if (length < minLineLength)
            {
                minLineLength = length;
            }
            if (length > maxLineLength)
            {
                maxLineLength = length;
            }
        }
    }
    //log_info << "minLineLength: " << minLineLength << " " << faceList.size();

    lengthSmallValue = smallValue * minLineLength;
    lengthVerySmallValue = verySmallValue * minLineLength;
    lengthLargeValue = largeValue * maxLineLength;
    lengthVeryLargeValue = veryLargeValue * maxLineLength;




    //============================== begin: creating serial mfem mesh ======================
    if (config->debugLevel > 1) log_info << "begin: creating mfem serial mesh";

    //mpi init at ParallelMPI.cpp
    //mfem::MPI_Session mfem::mpi;

    const char* device_config = "cpu";
    mfem::Device device(device_config);

    //create mfem vertices
    double* vertices;
    int num_vertices = pointList.size();
    vertices = new double[num_vertices * 3];


    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        auto& point = pointList[iPoint];
        vertices[3 * iPoint + 0] = point.position(0);
        vertices[3 * iPoint + 1] = point.position(1);
        vertices[3 * iPoint + 2] = point.position(2);
    }


    //create mfem element, for 2d, element is triangle
    std::vector<int> element_indices;
    std::vector<int> element_attributes;
    mfem::Geometry::Type element_type = mfem::Geometry::TRIANGLE;
    int num_elements = faceList.size();

    element_indices.resize(faceList.size() * 3);
    element_attributes.resize(faceList.size());
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        element_attributes[iFace] = 0;
        element_indices[iFace * 3 + 0] = faceList[iFace].pointList[0];
        element_indices[iFace * 3 + 1] = faceList[iFace].pointList[1];
        element_indices[iFace * 3 + 2] = faceList[iFace].pointList[2];

        element_attributes[iFace] = 1;
    }
    //** attribute for element
    int index = 10;
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        faceGroup.attributeIndex = index;
        index = index + 1;
        for (int i = 0; i < faceGroup.faceList.size(); i++)
        {
            //cout << "volumeGroup.volumeList[i]======="<<i<<"======" << volumeGroup.volumeList[i] << endl;
            element_attributes[faceGroup.faceList[i]] = faceGroup.attributeIndex;
        }
    }


    //create mfem boundary, for 2d, boudary element is segment
    std::vector<int> boundary_indices;
    std::vector<int> boundary_attributes;
    int num_boundary_elements = 0;
    mfem::Geometry::Type boundary_type = mfem::Geometry::SEGMENT;

    int attributeIndex = 1;
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto lineGroupName = iter->first;
        auto& lineGroup = iter->second;
        lineGroup.attributeIndex = attributeIndex;
        attributeIndex++;

        if (lineGroup.isInternal)
        {
            continue;
        }

        auto& lineElementList = lineGroup.lineElementList;
        for (int iLineElement = 0; iLineElement < lineElementList.size(); iLineElement++)
        {
            Element* lineElement = lineElementList[iLineElement];
            boundary_indices.push_back(lineElement->pointList[0]);
            boundary_indices.push_back(lineElement->pointList[1]);

            boundary_attributes.push_back(lineGroup.attributeIndex);
        }
        
    }
    num_boundary_elements = boundary_attributes.size();


    int dimension = config->dimension;
    int order = config->order;


   

    //============================== end: creating serial mfem mesh ======================



    // calculate lineList from mfem::Mesh
    //!!! mfem 2d 情况，GetNE得到的是面的个数，而GetNFaces是得不到面的个数的
    int numLine = mfemMesh->GetNEdges();
    if (config->debugLevel > 1) log_info << "serial mesh  vertex number: " << mfemMesh->GetNV();
    if (config->debugLevel > 1) log_info << "serial mesh  line number: " << mfemMesh->GetNEdges();
    if (config->debugLevel > 1) log_info << "serial mesh  face number: " << mfemMesh->GetNFaces();
    if (config->debugLevel > 1) log_info << "serial mesh  element number: " << mfemMesh->GetNE();
    lineList.resize(numLine);
    for (int iLine = 0; iLine < numLine; iLine++)
    {
        mfem::Array<int> vert;
        mfemMesh->GetEdgeVertices(iLine, vert);
        lineList[iLine].pointList[0] = vert[0];
        lineList[iLine].pointList[1] = vert[1];
    }

    numberPoint = pointList.size();
    numberLine = lineList.size();
    numberFace = faceList.size();
    numberVolume = volumeList.size();

    if (config->dimension == 3)
    {
        numberCell = numberVolume;
    }
    else if (config->dimension == 2)
    {
        numberCell = numberFace;
    }
    else
    {
        log_error << "dimension no right: " << config->dimension;
        exit(0);
    }

    
}


void Mesh::readGmshMesh()
{
    //log_info<<"11111";
    //gmsh::GmshMesh gmsh_Mesh(config->mesh_param.meshName);
    //log_info<<"22222";

    ////read point
    //pointList.resize(gmsh_Mesh.node_list.size());
    //for(int i = 0; i < pointList.size(); i++)
    //{
    //    pointList[i].position(0) = gmsh_Mesh.node_list[i].x;
    //    pointList[i].position(1) = gmsh_Mesh.node_list[i].y;
    //    pointList[i].position(2) = gmsh_Mesh.node_list[i].z;
    //}

    ///*
    //log_info<<"--------------------------------------------"<<endl;
    //for(int i = 0; i < vertexes.size(); i++)
    //{
    //    log_info<<"vertex "<<i<<" "<<vertexes[i].x[0]<<" "<<vertexes[i].x[1]<<endl;
    //}
    //*/

    //log_info<<"33333  point 10: "<<pointList[10].position(0)<<" "<<pointList[10].position(1) <<" "<<pointList[10].position(2);
    ////read pointGroup, lineGroup, faceGroup, volumeGroup
    //int physicalGroupNumber[4];
    //for(auto it = gmsh_Mesh.physical_group_map.begin(); it != gmsh_Mesh.physical_group_map.end(); it++)
    //{
    //    auto& physical_group = it->second;
    //    if(physical_group.dim == 0)
    //    {
    //        PointGroup pointGroup;
    //        pointGroup.name = physical_group.name;
    //        pointGroup.pointList.resize(0);
    //        pointGroupMap[pointGroup.name] = pointGroup;
    //    }
    //    else if(physical_group.dim == 1)
    //    {
    //        LineGroup lineGroup;
    //        lineGroup.name = physical_group.name;
    //        lineGroup.lineList.resize(0);
    //        lineGroupMap[lineGroup.name] = lineGroup;
    //    }
    //    else if(physical_group.dim == 2)
    //    {
    //        FaceGroup faceGroup;
    //        faceGroup.name = physical_group.name;
    //        faceGroup.faceList.resize(0);
    //        faceGroupMap[faceGroup.name] = faceGroup;
    //    }
    //    else if(physical_group.dim == 3)
    //    {
    //        VolumeGroup volumeGroup;
    //        volumeGroup.name = physical_group.name;
    //        volumeGroup.volumeList.resize(0);
    //        volumeGroupMap[volumeGroup.name] = volumeGroup;
    //    }
    //}


    //log_info<<"44444 volumeList size: " << gmsh_Mesh.element_list.size();
    ////read elements, including: point, line, face, volume
    ////read physical group togather
    //for(int i = 0; i < gmsh_Mesh.element_list.size(); i++)
    //{
    //    auto& element = gmsh_Mesh.element_list[i];
    //    //log_info<<"44444 element number " << i << " element type: " << element.element_type;

    //    if(element.element_type == 1)
    //    {
    //        Line line;
    //        //gmsh index start from 1
    //        line.pointList[0] = element.node_tag_list[0] - 1;
    //        line.pointList[1] = element.node_tag_list[1] - 1;
    //        lineList.push_back(line);

    //        int entity_tag = element.entity_tag;
    //        auto& physical_tag_list = gmsh_Mesh.curve_entity_map[entity_tag].physical_tag_list;
    //        for(int j = 0; j < physical_tag_list.size(); j++)
    //        {
    //            int physical_tag = physical_tag_list[j];
    //            std::string name = gmsh_Mesh.physical_group_map[physical_tag].name;
    //            LineGroup& lineGroup = lineGroupMap[name];
    //            lineGroup.lineList.push_back(lineList.size()-1);
    //        }
    //    }
    //    else if(element.element_type == 2)
    //    {
    //        //log_info<<"----------- face list size " << faceList.size() << " node_tag_list size: " << element.node_tag_list.size();
    //        Face face;

    //        //gmsh index start from 1, 3 point direction is wiseclock
    //        //easypic3d Mesh trianle, 3 point direction is wiseclock
    //        face.pointList[0] = element.node_tag_list[0] - 1;
    //        face.pointList[1] = element.node_tag_list[2] - 1;
    //        face.pointList[2] = element.node_tag_list[1] - 1;

    //        
    //        faceList.push_back(face);

    //        int entity_tag = element.entity_tag;
    //        //log_info<<"----------- entity_tag " << entity_tag;

    //        auto& physical_tag_list = gmsh_Mesh.surface_entity_map[entity_tag].physical_tag_list;
    //        for(int j = 0; j < physical_tag_list.size(); j++)
    //        {
    //            //log_info<<"----------- physical name j " << j;
    //            int physical_tag = physical_tag_list[j];
    //            std::string name = gmsh_Mesh.physical_group_map[physical_tag].name;

    //            //log_info<<"----------- physical name " << name;
    //            FaceGroup& faceGroup = faceGroupMap[name];
    //            faceGroup.faceList.push_back(faceList.size()-1);
    //        }
    //    }
    //    else if(element.element_type == 4)
    //    {
    //        Volume volume;

    //        //gmsh index start from 1, 3 point direction is wiseclock
    //        //easypic3d Mesh trianle, 3 point direction is wiseclock
    //        volume.pointList[0] = element.node_tag_list[0] - 1;
    //        volume.pointList[1] = element.node_tag_list[3] - 1;
    //        volume.pointList[2] = element.node_tag_list[2] - 1;
    //        volume.pointList[3] = element.node_tag_list[1] - 1;
    //        volumeList.push_back(volume);


    //        int entity_tag = element.entity_tag;
    //        auto& physical_tag_list = gmsh_Mesh.curve_entity_map[entity_tag].physical_tag_list;
    //        for(int j = 0; j < physical_tag_list.size(); j++)
    //        {
    //            int physical_tag = physical_tag_list[j];
    //            std::string name = gmsh_Mesh.physical_group_map[physical_tag].name;
    //            VolumeGroup& volumeGroup = volumeGroupMap[name];
    //            volumeGroup.volumeList.push_back(volumeList.size()-1);
    //        }
    //    }
    //}
    //log_info << "[Mesh.cpp::readGmshMesh] reading line done";



}


void Mesh::setMaterial()
{
    for (auto iter = volumeGroupMap.begin(); iter != volumeGroupMap.end(); iter++)
    {
        auto& volumeGroup = iter->second;

        //get volumeGroup material
        auto volumeGroupToMaterial = config->mesh_param.physicalGroupToMaterialMap.find(volumeGroup.name);
        if (volumeGroupToMaterial == config->mesh_param.physicalGroupToMaterialMap.end())
        {
            //volumeGroupName没有对应的material,则设置材料为真空
            volumeGroup.material.relativePermittivity = 1.0;
        }
        else
        {
            std::string materialName = volumeGroupToMaterial->second;
            auto material = config->materialMap.find(materialName);
            if (material == config->materialMap.end())
            {
                //没有对应的material
                log_error << "no such material in database: " << materialName;
            }
            else
            {
                volumeGroup.material = material->second;
            }

        }
    }
}


void Mesh::setMaterial2d()
{
    
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto& faceGroup = iter->second;

        //get faceGroup material
        auto physicalGroupToMaterial = config->mesh_param.physicalGroupToMaterialMap.find(faceGroup.name);
        if (physicalGroupToMaterial == config->mesh_param.physicalGroupToMaterialMap.end())
        {
            //physicalGroupName没有对应的material,则设置材料为真空
            faceGroup.material.relativePermittivity = 1.0;
        }
        else
        {
            std::string materialName = physicalGroupToMaterial->second;
            auto material = config->materialMap.find(materialName);
            if (material == config->materialMap.end())
            {
                //没有对应的material
                log_error << "no such material in database: " << materialName;
            }
            else
            {
                faceGroup.material = material->second;
            }

        }
    }
    
}



//cgns和gmsh的单元是一起存储的，包括line face volume，需要将line face volume分开存储，重新编号，因为EasyPIC中faceGroup中是face的编号
//而且cgns和gmsh存储的face只是边界face，没有内部face，还需要得到内部和边界所有face，然后将faceGroup中边界face编号转换成全局face编号
//cal connectivity
void Mesh::calConnectivity()
{
    if (config->debugLevel > 1) log_info << "begin: calConnectivity";
    //faceMap作为临时变量，用完之后内存会自动释放，避免占用过多内存
    std::map<Face, int> faceMap;
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        faceMap[faceList[iFace]] = iFace;
    }

    //log_info << "88888";
    //不再需要额外确定volumeGroup中单元在volumeList中的序号，因为在readMesh的时候，单元列表里包含了所有体单元，而不像面单元那样，
    //gmsh cgns等格式网格只存储了边界面单元。因此faceGroup需要重新计算序号，而体单元不需要重新计算。

    //calculte volumeToFaceList and faceToVolumeList
    volumeToFaceList.resize(volumeList.size());
    faceToVolumeList.resize(faceList.size());
    for (int iFace = 0; iFace < faceToVolumeList.size(); iFace++)
    {
        FaceToVolume& faceToVolume = faceToVolumeList[iFace];
        //for boundary face, has only one volume, the other is set to -1
        faceToVolume.volumeList[0] = -1;
        faceToVolume.volumeList[1] = -1;
    }
    for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        const Volume& volume = volumeList[iVolume];
        VolumeToFace& volumeToFace = volumeToFaceList[iVolume];
        for (int iFace = 0; iFace < 4; iFace++)
        {
            Face faceLocal;
            for (int iPoint = 0; iPoint < 3; iPoint++)
            {
                //点按顺时针存取，faceLocal按逆序存储
                //faceLocal.pointList[iPoint] = volume.pointList[(iFace+iPoint+1)%4] //shun xu
                faceLocal.pointList[iPoint] = volume.pointList[(iFace - iPoint + 3) % 4]; //ni xu
            }
            auto iter = findFace(faceMap, faceLocal);
            if (iter != faceMap.end())
            {
                //has found
                int iFaceGlobal = iter->second;
                volumeToFace.faceList[iFace] = iFaceGlobal;
                if (faceToVolumeList[iFaceGlobal].volumeList[0] == -1)
                {
                    faceToVolumeList[iFaceGlobal].volumeList[0] = iVolume;
                }
                else
                {
                    faceToVolumeList[iFaceGlobal].volumeList[1] = iVolume;
                }
                
            }
            else
            {
                //not found
                log_error << "faceToVolume cal error";
            }

        }

    }

    //calculate volumeToFaceToVolume from ListvolumeToFaceList and faceToVolumeList
    volumeToFaceToVolumeList.resize(volumeList.size());
    for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        VolumeToFace& volumeToFace = volumeToFaceList[iVolume];
        VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];
        for (int iFace = 0; iFace < 4; iFace++)
        {
            volumeToFaceToVolume.faceList[iFace] = volumeToFace.faceList[iFace];
        }

        for (int iFace = 0; iFace < 4; iFace++)
        {
            int iFaceGlobal = volumeToFaceToVolume.faceList[iFace];
            FaceToVolume& faceToVolume = faceToVolumeList[iFaceGlobal];
            if (faceToVolume.volumeList[0] == iVolume)
            {
                if (faceToVolume.volumeList[1] == -1)
                {
                    //this face is boundary
                }
                else
                {
                    volumeToFaceToVolume.volumeList[iFace] = faceToVolume.volumeList[1];
                }

            }
            else
            {
                volumeToFaceToVolume.volumeList[iFace] = faceToVolume.volumeList[0];
            }

        }

    }

    //cal VolumeToFaceToVolume other info: faceCenterPoint, normalList
    if (config->debugLevel > 1) log_info << "begin: cal VolumeToFaceToVolume other info";
    for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        Volume& volume = volumeList[iVolume];
        Point& point0 = pointList[volume.pointList[0]];
        Point& point1 = pointList[volume.pointList[1]];
        Point& point2 = pointList[volume.pointList[2]];
        Point& point3 = pointList[volume.pointList[3]];


        //cal volumeToFaceToVolume::centerPoint
        VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];
        for (int i = 0; i < 3; i++)
        {
            volumeToFaceToVolume.centerPoint.position(i) = 0.0;
        }
        for (int i = 0; i < 3; i++)
        {
            volumeToFaceToVolume.centerPoint.position(i) += point0.position(i) * 0.25;
            volumeToFaceToVolume.centerPoint.position(i) += point1.position(i) * 0.25;
            volumeToFaceToVolume.centerPoint.position(i) += point2.position(i) * 0.25;
            volumeToFaceToVolume.centerPoint.position(i) += point3.position(i) * 0.25;
        }

        //cal volumeToFaceToVolume::faceCenterPointList
        for (int iFace = 0; iFace < 4; iFace++)
        {
            Point& faceCenterPoint = volumeToFaceToVolume.faceCenterPointList[iFace];
            Face& face = faceList[volumeToFaceToVolume.faceList[iFace]];
            for (int i = 0; i < 3; i++)
            {
                faceCenterPoint.position(i) = 0.0;
            }
            for (int i = 0; i < 3; i++)
            {
                for (int iPoint = 0; iPoint < 3; iPoint++)
                {
                    Point& point = pointList[face.pointList[iPoint]];
                    faceCenterPoint.position(i) += point.position(i) / 3.0;
                }
            }
        }

        //cal volumeToFaceToVolume::normalList(this face out normal)
        /*for (int iFace = 0; iFace < 4; iFace++)
        {
            Vector3d& normal = volumeToFaceToVolume.normalList[iFace];
            Face& face = faceList[volumeToFaceToVolume.faceList[iFace]];

            Point& point0 = pointList[face.pointList[0]];
            Point& point1 = pointList[face.pointList[1]];
            Point& point2 = pointList[face.pointList[2]];

            Vector3d A;
            Vector3d B;
            Vector3d C;

            Vector3d ab;
            Vector3d ac;
            Vector3d abcNormal;
            Vector3d volumeCenterToFaceCenter;

            volumeCenterToFaceCenter(0) = volumeToFaceToVolume.faceCenterPointList[iFace].position(0) - volumeToFaceToVolume.centerPoint.position(0);
            volumeCenterToFaceCenter(1) = volumeToFaceToVolume.faceCenterPointList[iFace].position(1) - volumeToFaceToVolume.centerPoint.position(1);
            volumeCenterToFaceCenter(2) = volumeToFaceToVolume.faceCenterPointList[iFace].position(2) - volumeToFaceToVolume.centerPoint.position(2);


            for (int i = 0; i < 3; i++)
            {
                A(i) = point0.position(i);
                B(i) = point1.position(i);
                C(i) = point2.position(i);
            }
            ab = B - A;
            ac = C - A;
            abcNormal = ab.cross(ac);
            abcNormal.normalize();

            if (volumeCenterToFaceCenter.dot(abcNormal) > 0.0)
            {
                normal(0) = abcNormal[0];
                normal(1) = abcNormal[1];
                normal(2) = abcNormal[2];
            }
            else
            {
                normal(0) = -abcNormal[0];
                normal(1) = -abcNormal[1];
                normal(2) = -abcNormal[2];
            }



        }*/

    }
    //log_info << "end: cal VolumeToFaceToVolume other info";


    //set faceCenterPoint for face in faceList
    //log_info << "begin: set faceCenterPoint for face in faceList";
    for (int iVolume = 0; iVolume < volumeToFaceToVolumeList.size(); iVolume++)
    {
        auto& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];
        for (int iFace = 0; iFace < 4; iFace++)
        {
            int iFaceGlobal = volumeToFaceToVolume.faceList[iFace];
            faceList[iFaceGlobal].faceCenterPoint = volumeToFaceToVolume.faceCenterPointList[iFace];
        }
    }


    //==================================== begin: calculate for faceGroup =================================
    //确定faceGroup中单元在faceList中的序号
    //gmsh的只存了边界的face，faceGroup里边界的face序号也是边界face的，
    //EasyPIC里存的是所有face（包括边界和内部），下面是将faceGroup里face序号转换到全局face序号
    //face中点的顺序有六种可能，因此需要六次判断
    //log_info << "begin: set face index of faceGroup";
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        std::string faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        faceGroup.faceList.resize(faceGroup.faceElementList.size());
        for (int i = 0; i < faceGroup.faceList.size(); i++)
        {
            int iFaceLocal;
            Face face;
            face.pointList[0] = faceGroup.faceElementList[i]->pointList[0];
            face.pointList[1] = faceGroup.faceElementList[i]->pointList[1];
            face.pointList[2] = faceGroup.faceElementList[i]->pointList[2];

            auto iter = findFace(faceMap, face);
            if (iter != faceMap.end())
            {
                //has found
                int iFace = iter->second;
                faceGroup.faceList[i] = iFace;

                //log_info << "global faceGroup------------" << faceGroupName << " " << i << " " << iFace;
            }
            else
            {
                log_fatal << "not found, faceGroup: " << faceGroupName;
            }

        }
    }


    //calculate faceGroup of volumeGroup's surface
    //log_info << "begin: calculate faceGroup of volumeGroup's surface";
    //!!!this process should be executed in geometry module of GUI
    //for (auto iter = volumeGroupMap.begin(); iter != volumeGroupMap.end(); iter++)
    //{
    //    string faceGroupName = iter->first + "Surface";
    //    auto& volumeGroup = iter->second;

    //    vector<int> faceListTemp;

    //    map<int, int> volumeMapTemp;
    //    for (int i = 0; i < volumeGroup.volumeList.size(); i++)
    //    {
    //        int iVolume = volumeGroup.volumeList[i];
    //        volumeMapTemp[iVolume] = iVolume;
    //    }

    //    for (int i = 0; i < volumeGroup.volumeList.size(); i++)
    //    {
    //        int iVolume = volumeGroup.volumeList[i];
    //        VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];

    //        for (int j = 0; j < 4; j++)
    //        {
    //            int iFace = volumeToFaceToVolume.faceList[j];
    //            FaceToVolume& faceToVolume = faceToVolumeList[iFace];
    //            if (faceToVolume.volumeList[1] == -1)
    //            {
    //                continue;
    //            }

    //            int numIn = 0;
    //            int numWhich = -1;
    //            auto iter = volumeMapTemp.find(faceToVolume.volumeList[0]);
    //            if (iter != volumeMapTemp.end())
    //            {
    //                numIn++;
    //                numWhich = 0;
    //            }
    //            iter = volumeMapTemp.find(faceToVolume.volumeList[1]);
    //            if (iter != volumeMapTemp.end())
    //            {
    //                numIn++;
    //                numWhich = 1;
    //            }

    //            if (numIn == 1)
    //            {
    //                faceListTemp.push_back(iFace);
    //                faceList[iFace].normal = volumeToFaceToVolume.normalList[j];
    //                faceList[iFace].isInGroup = true;

    //                //faceToVolume.volumeList[]的第一个volume属于vacuum region
    //                if (numWhich == 0)
    //                {
    //                    int iVolumeTemp = faceToVolume.volumeList[0];
    //                    faceToVolume.volumeList[0] = faceToVolume.volumeList[1];
    //                    faceToVolume.volumeList[1] = iVolumeTemp;
    //                }
    //            }
    //        }
    //    }

    //    FaceGroup faceGroup;
    //    faceGroup.name = faceGroupName;
    //    faceGroup.faceList = faceListTemp;
    //    faceGroupMap[faceGroupName] = faceGroup;
    //}


    //reorder faceToVolume.volumeList for objectSurface (the faceGroup corresponding to volumeGroup)
    //make sure faceToVolume.volumeList[0] is in vacuum region
    for (auto iter = volumeGroupMap.begin(); iter != volumeGroupMap.end(); iter++)
    {
        std::string faceGroupName = iter->first + "Surface";
        auto& volumeGroup = iter->second;

        std::map<int, int> volumeMapTemp;
        for (int i = 0; i < volumeGroup.volumeList.size(); i++)
        {
            int iVolume = volumeGroup.volumeList[i];
            volumeMapTemp[iVolume] = iVolume;
        }

        for (int i = 0; i < volumeGroup.volumeList.size(); i++)
        {
            int iVolume = volumeGroup.volumeList[i];
            VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];

            for (int j = 0; j < 4; j++)
            {
                int iFace = volumeToFaceToVolume.faceList[j];
                FaceToVolume& faceToVolume = faceToVolumeList[iFace];
                if (faceToVolume.volumeList[1] == -1)
                {
                    continue;
                }

                int numIn = 0;
                int numWhich = -1;
                auto iter = volumeMapTemp.find(faceToVolume.volumeList[0]);
                if (iter != volumeMapTemp.end())
                {
                    numIn++;
                    numWhich = 0;
                }
                iter = volumeMapTemp.find(faceToVolume.volumeList[1]);
                if (iter != volumeMapTemp.end())
                {
                    numIn++;
                    numWhich = 1;
                }

                if (numIn == 1)
                {
                    faceList[iFace].normal = volumeToFaceToVolume.normalList[j];
                    faceList[iFace].isInGroup = true;

                    //faceToVolume.volumeList[]的第一个volume属于vacuum region
                    if (numWhich == 0)
                    {
                        int iVolumeTemp = faceToVolume.volumeList[0];
                        faceToVolume.volumeList[0] = faceToVolume.volumeList[1];
                        faceToVolume.volumeList[1] = iVolumeTemp;
                    }
                }
            }
        }
    }


    //calculate pointList of faceGroup
    //log_info << "begin: calculate pointList of faceGroup";
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto& faceGroupName = iter->first;
        auto& faceGroup = iter->second;

        std::set<int> pointSet;
        for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
        {
            auto& face = faceList[faceGroup.faceList[iFace]];
            for (int iPoint = 0; iPoint < 3; iPoint++)
            {
                pointSet.insert(face.pointList[iPoint]);
            }
        }

        faceGroup.pointList.clear();
        for (auto iterPoint = pointSet.begin(); iterPoint != pointSet.end(); iterPoint++)
        {
            faceGroup.pointList.push_back(*iterPoint);
        }
    }

    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto& faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        faceGroup.isPeriodic = false;

        faceGroup.numberFace = faceGroup.faceList.size();
        faceGroup.numberFaceSerial = faceGroup.numberFace;
    }



    //check point number for periodic boundary faceGroup
    //reorder periodic faceGroup pointList
    //log_info << "begin: reorder periodic faceGroup pointList";
    for (int iFaceGroup = 0; iFaceGroup < config->mesh_param.mesh_periodic_edge0.size(); iFaceGroup++)
    {
        std::string faceName0 = config->mesh_param.mesh_periodic_edge0[iFaceGroup];
        std::string faceName1 = config->mesh_param.mesh_periodic_edge1[iFaceGroup];
        auto& faceGroup0 = faceGroupMap.find(faceName0)->second;
        auto& faceGroup1 = faceGroupMap.find(faceName1)->second;

        //log_info << "faceName0 and faceName1: " << faceName0 << " " << faceGroup0.pointList.size() << " " << faceName1 << " " << faceGroup1.pointList.size();
        if (faceGroup0.pointList.size() != faceGroup1.pointList.size())
        {
            log_error << "periodic boundary faceGroup do not match: " << faceName0 << " " << faceGroup0.pointList.size() << " " << faceName1 << " " << faceGroup1.pointList.size();
            exit(0);
        }

        std::vector<int> pointListNew1;
        pointListNew1.resize(faceGroup0.pointList.size());
        for (int i = 0; i < faceGroup0.pointList.size(); i++)
        {
            int iPoint0 = faceGroup0.pointList[i];
            Point& point0 = pointList[iPoint0];

            double minLength = maxGeometryLength;
            int iPoint;
            for (int j = 0; j < faceGroup1.pointList.size(); j++)
            {
                int iPoint1 = faceGroup1.pointList[j];
                Point& point1 = pointList[iPoint1];
                double lineLength = calLineLength(point0, point1);
                if (lineLength < minLength)
                {
                    minLength = lineLength;
                    iPoint = iPoint1;
                }
            }
            pointListNew1[i] = iPoint;
            //log_info << "iPoint: " << iPoint;

        }
        faceGroup1.pointList = pointListNew1;

        //log_info << "pointListNew1.size: " << pointListNew1.size();
        //log_info << "faceName0 and faceName1: " << faceName0 << " " << faceGroup0.pointList.size() << " " << faceName1 << " " << faceGroup1.pointList.size();

        //log_info << "2222";
        //calculate periodicTranslationVector
        auto& point0 = pointList[faceGroup0.pointList[0]];
        auto& point1 = pointList[faceGroup1.pointList[0]];
        //log_info << "position: " << point1.position(0);

        //log_info << "2222 aaa";
        faceGroup0.periodicTranslationVector(0) = point1.position(0) - point0.position(0);
        //log_info << "2222 bbb";
        faceGroup0.periodicTranslationVector(1) = point1.position(1) - point0.position(1);
        faceGroup0.periodicTranslationVector(2) = point1.position(2) - point0.position(2);
        faceGroup1.periodicTranslationVector(0) = -faceGroup0.periodicTranslationVector(0);
        faceGroup1.periodicTranslationVector(1) = -faceGroup0.periodicTranslationVector(1);
        faceGroup1.periodicTranslationVector(2) = -faceGroup0.periodicTranslationVector(2);
        //log_info << "3333";
        faceGroup0.isPeriodic = true;
        faceGroup1.isPeriodic = true;

    }




    //reorder periodic faceGroup faceList
    //log_info << "begin: reorder periodic faceGroup faceList";
    for (int iFaceGroup = 0; iFaceGroup < config->mesh_param.mesh_periodic_edge0.size(); iFaceGroup++)
    {
        std::string faceName0 = config->mesh_param.mesh_periodic_edge0[iFaceGroup];
        std::string faceName1 = config->mesh_param.mesh_periodic_edge1[iFaceGroup];
        auto& faceGroup0 = faceGroupMap.find(faceName0)->second;
        auto& faceGroup1 = faceGroupMap.find(faceName1)->second;
        //log_info << "++++++++++++++++++++++++++++++++++++++++++++++++++++++";
        //log_info << faceName0 << " " << faceName1;


        std::vector<int> faceListNew1;
        faceListNew1.resize(faceGroup0.faceList.size());

        for (int i = 0; i < faceGroup0.faceList.size(); i++)
        {
            int iFace0 = faceGroup0.faceList[i];
            Point& point0 = faceList[iFace0].faceCenterPoint;


            double minLength = maxGeometryLength;
            int iFace;
            for (int j = 0; j < faceGroup1.faceList.size(); j++)
            {
                int iFace1 = faceGroup1.faceList[j];
                Point& point1 = faceList[iFace1].faceCenterPoint;
                double lineLength = calLineLength(point0, point1);
                if (lineLength < minLength)
                {
                    minLength = lineLength;
                    iFace = iFace1;
                }
            }
            faceListNew1[i] = iFace;

            faceList[iFace0].periodicFaceIndex = iFace;
            faceList[iFace].periodicFaceIndex = iFace0;

        }
        faceGroup1.faceList = faceListNew1;
    }
    //log_info << "end: reorder periodic faceGroup faceList";




    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        faceList[iFace].isInGroup = false;
    }

    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        //log_info << faceGroupName << " face number: " << faceGroup.faceList.size();
        for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
        {
            //log_info << " ----------- face id: " << faceGroup.faceList[iFace];
            faceList[faceGroup.faceList[iFace]].isInGroup = true;
            faceList[faceGroup.faceList[iFace]].faceGroupName = faceGroupName;
        }

    }


    numBoundaryFace = 0;
    for(int iFace = 0; iFace < faceToVolumeList.size(); iFace++)
    {
        auto& faceToVolume = faceToVolumeList[iFace];
        if(faceToVolume.volumeList[1] == -1)
        {
            numBoundaryFace++;
        }
    }
    if (config->debugLevel > 1) log_info << "end: calculate for faceGroup" <<std::endl;
    //==================================== end: calculate for faceGroup =================================







    //cal area of face
    //两个边的向量叉乘除以2
    //log_info << "cal area of face";
    //for(int iFace = 0; iFace < faceList.size(); iFace++)
    //{
    //    Face& face = faceList[iFace];
    //    Point& point0 = pointList[face.pointList[0]];
    //    Point& point1 = pointList[face.pointList[1]];
    //    Point& point2 = pointList[face.pointList[2]];

    //    Eigen::Vector3d vector0;
    //    vector0(0) = point1.position(0) - point0.position(0);
    //    vector0(1) = point1.position(1) - point0.position(1);
    //    vector0(2) = point1.position(2) - point0.position(2);

    //    Eigen::Vector3d vector1;
    //    vector1(0) = point2.position(0) - point0.position(0);
    //    vector1(1) = point2.position(1) - point0.position(1);
    //    vector1(2) = point2.position(2) - point0.position(2);

    //    face.area = vector0.cross(vector1).norm() / 2.0;
    //    //log_info << "face area: " << iFace << " " << face.area<< " " << face.pointList[0] << " " << face.pointList[1] << " " << face.pointList[2];
    //}


    //cal volume of Volume(tetrahedron)
    //log_info << "start calculate volume of tetrahedron";
    minVol = 1.0e99;
    maxVol = 0.0;
    totalVol = 0.0;
    /*for(int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        Volume& volume = volumeList[iVolume];
        Point& point0 = pointList[volume.pointList[0]];
        Point& point1 = pointList[volume.pointList[1]];
        Point& point2 = pointList[volume.pointList[2]];
        Point& point3 = pointList[volume.pointList[3]];

        Eigen::Matrix3d matrix;
        matrix(0, 0) = point1.position(0) - point0.position(0);
        matrix(0, 1) = point2.position(0) - point0.position(0);
        matrix(0, 2) = point3.position(0) - point0.position(0);
        matrix(1, 0) = point1.position(1) - point0.position(1);
        matrix(1, 1) = point2.position(1) - point0.position(1);
        matrix(1, 2) = point3.position(1) - point0.position(1);
        matrix(2, 0) = point1.position(2) - point0.position(2);
        matrix(2, 1) = point2.position(2) - point0.position(2);
        matrix(2, 2) = point3.position(2) - point0.position(2);

        volume.vol = abs(matrix.determinant()) / 6.0;

        totalVol += volume.vol;

        if(volume.vol > maxVol)
        {
            maxVol = volume.vol;
        }

        if(minVol == 0.0)
        {
            minVol = volume.vol;
        }
        if(volume.vol < minVol)
        {
            minVol = volume.vol;
        }

    }*/
    averageVol = totalVol / volumeList.size();





    //cal dualVolumeVolList
    //ref: 2018_a parallel electrostatic particle in cell method on unstructured tetrahedral grids for large scale bounded collisionless plasma simulatins
    //每个体中对应的点有一个六面体，计算这个六面体的面积加到这个点对应的dualVolume上，计算方法：六面体六个面分别和六面体中心构成六个四棱锥，每个四棱锥分成两个四面体，计算四面体体积然后相加
    //log_info << "start calculate dual volume";
    dualVolumeVolList.resize(pointList.size());
    for(int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        dualVolumeVolList[iPoint] = 0.0;
    }
    for(int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    {
        Volume& volume = volumeList[iVolume];
        for(int iPoint = 0; iPoint < 4; iPoint++)
        {
            int iPointGlobal;
            loc2glb(iVolume, iPoint, iPointGlobal);

            Point& point0 = pointList[volume.pointList[(iPoint+0)%4]];
            Point& point1 = pointList[volume.pointList[(iPoint+1)%4]];
            Point& point2 = pointList[volume.pointList[(iPoint+2)%4]];
            Point& point3 = pointList[volume.pointList[(iPoint+3)%4]];

            Point& pointVolumeCenter = volumeToFaceToVolumeList[iVolume].centerPoint;

            Point pointFaceCenter0;
            Point pointFaceCenter1;
            Point pointFaceCenter2;
            for(int i = 0; i < 3; i++)
            {
                pointFaceCenter0.position(i) = (point0.position(i) + point1.position(i) + point2.position(i)) / 3.0;
                pointFaceCenter1.position(i) = (point0.position(i) + point2.position(i) + point3.position(i)) / 3.0;
                pointFaceCenter2.position(i) = (point0.position(i) + point3.position(i) + point1.position(i)) / 3.0;
            }

            Point pointLineCenter0;
            Point pointLineCenter1;
            Point pointLineCenter2;
            for(int i = 0; i < 3; i++)
            {
                pointLineCenter0.position(i) = (point0.position(i) + point1.position(i)) / 2.0;
                pointLineCenter1.position(i) = (point0.position(i) + point2.position(i)) / 2.0;
                pointLineCenter2.position(i) = (point0.position(i) + point3.position(i)) / 2.0;
            }

            //hex point order: point0, l0, f0, l1, l2, f2, volumeCenter, f1
            Point hexCenterPoint;
            for(int i = 0; i < 3; i++)
            {
                hexCenterPoint.position(i) = (point0.position(i) + pointVolumeCenter.position(i) + pointFaceCenter0.position(i) + pointFaceCenter1.position(i) + pointFaceCenter2.position(i) 
                                           + pointLineCenter0.position(i) + pointLineCenter1.position(i) + pointLineCenter2.position(i)) / 8.0;
            }
            

            //cal face: point0,l0, f0, l1
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter0, pointLineCenter1, point0);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter0, pointLineCenter1, pointFaceCenter0);

            //cal face: l2, f2, volumeCenter, f1
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointFaceCenter1, pointFaceCenter2, pointLineCenter2);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointFaceCenter1, pointFaceCenter2, pointVolumeCenter);

            //cal face: point0,l0, f2, l2
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter0, pointLineCenter2, point0);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter0, pointLineCenter2, pointFaceCenter2);

            //cal face: l0, f0, volumeCenter, f2
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointFaceCenter0, pointFaceCenter2, pointLineCenter0);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointFaceCenter0, pointFaceCenter2, pointVolumeCenter);

            //cal face: f0, l1, f1, volumeCenter
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter1, pointVolumeCenter, pointFaceCenter0);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, pointLineCenter1, pointVolumeCenter, pointFaceCenter1);

            //cal face: l1, point0, l2, f1
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, point0, pointFaceCenter1, pointLineCenter1);
            dualVolumeVolList[iPointGlobal] += calTetVolume(hexCenterPoint, point0, pointFaceCenter1, pointLineCenter2);
        }
    }

    minDualVolumeVol = 1.0e99;
    maxDualVolumeVol = 0.0;

    for(int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        double vol = dualVolumeVolList[iPoint];
        if(vol > maxDualVolumeVol)
        {
            maxDualVolumeVol = vol;
        }
        if(vol < minDualVolumeVol)
        {
            minDualVolumeVol = vol;
        }
    }






    ////calculate barycentric coordinate matrix 
    ////log_info << "begin calculate barycentric coordinate matrix ";
    //barycentric_coordinate_A_inv.resize(volumeList.size());
    //TTensorList.resize(volumeList.size());
    //detAList.resize(volumeList.size());
    //for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    //{
    //    Volume& volume = volumeList[iVolume];

    //    //eigen bug: inverse of MatrixXd is OK, but inverse of Matrix4d is wrong, but data store type of barycentric_coordinate_A_inv use Matrix4d, which 
    //    //should be OK.
    //    Eigen::MatrixXd matrix(4, 4);
    //    //MatrixXd matrix1(4,4);
    //    //Matrix4d matrix;
    //    for (int iPoint = 0; iPoint < 4; iPoint++)
    //    {
    //        Point& point = pointList[volume.pointList[iPoint]];

    //        matrix(0, iPoint) = point.position(0);
    //        matrix(1, iPoint) = point.position(1);
    //        matrix(2, iPoint) = point.position(2);
    //        matrix(3, iPoint) = 1.0;
    //    }
    //    //matrix1 = matrix.inverse();
    //    barycentric_coordinate_A_inv[iVolume] = matrix.inverse();

    //    Eigen::Matrix<double, 4, 3> matrix43;
    //    Eigen::Vector3d ab, ac, ad, bc, bd, crossVector;

    //    Eigen::Map<Eigen::Vector3d> a(pointList[volume.pointList[0]].position);
    //    Eigen::Map<Eigen::Vector3d> b(pointList[volume.pointList[1]].position);
    //    Eigen::Map<Eigen::Vector3d> c(pointList[volume.pointList[2]].position);
    //    Eigen::Map<Eigen::Vector3d> d(pointList[volume.pointList[3]].position);

    //    ab = b - a;
    //    ac = c - a;
    //    ad = d - a;
    //    bc = c - b;
    //    bd = d - b;

    //    crossVector = bd.cross(bc);
    //    matrix43(0, 0) = crossVector(0);
    //    matrix43(0, 1) = crossVector(1);
    //    matrix43(0, 2) = crossVector(2);

    //    crossVector = ac.cross(ad);
    //    matrix43(1, 0) = crossVector(0);
    //    matrix43(1, 1) = crossVector(1);
    //    matrix43(1, 2) = crossVector(2);

    //    crossVector = ad.cross(ab);
    //    matrix43(2, 0) = crossVector(0);
    //    matrix43(2, 1) = crossVector(1);
    //    matrix43(2, 2) = crossVector(2);

    //    crossVector = ab.cross(ac);
    //    matrix43(3, 0) = crossVector(0);
    //    matrix43(3, 1) = crossVector(1);
    //    matrix43(3, 2) = crossVector(2);

    //    TTensorList[iVolume] = matrix43;

    //    detAList[iVolume] = ab.dot(ac.cross(ad));

    //}
    ////log_info << "end calculate barycentric coordinate matrix ";


     //计算lineGroup或者faceGroup的总面积
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto& faceGroupName = iter->first;
        auto& faceGroup = iter->second;

        double totalAreaTemp = 0.0;
        for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
        {
            auto& face = faceList[faceGroup.faceList[iFace]];
            totalAreaTemp += face.area;
        }

        faceGroup.totalArea = totalAreaTemp;
    }

    outputMeshInfo();

}



void Mesh::calConnectivity2d()
{

    //for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    //{
    //    auto& point = pointList[iPoint];
    //    log_info << "point position: " << iPoint << " " << point.position(0) << " " << point.position(1) << " " << point.position(2);
    //}

    //cal line.centPoint
    for (int iLine = 0; iLine < lineList.size(); iLine++)
    {
        auto& line = lineList[iLine];
        auto& point0 = pointList[line.pointList[0]];
        auto& point1 = pointList[line.pointList[1]];

        for (int i = 0; i < 3; i++)
        {
            line.lineCenterPoint.position(i) = 0.5 * (point0.position(i) + point1.position(i));
        }
        //log_info << "line center point position: " << iLine << " " << line.lineCenterPoint.position(0) << " " << line.lineCenterPoint.position(1) << " " << line.lineCenterPoint.position(2);
    }

    
    //lineMap作为临时变量，用完之后内存会自动释放，避免占用过多内存
    if (config->debugLevel > 1) log_info << "begin creating lineMap";
    std::map<Line, int> lineMap;
    for (int iLine = 0; iLine < lineList.size(); iLine++)
    {
        lineMap[lineList[iLine]] = iLine;
    }


    //gmsh cgns等格式网格只存储了边界线单元。因此lineGroup需要重新计算序号，而面单元不需要重新计算。

    //calculte volumeToFaceList and faceToVolumeList
    if (config->debugLevel > 1) log_info << "begin: creating faceTolineList and lineToFaceList";
    faceToLineList.resize(faceList.size());
    lineToFaceList.resize(lineList.size());
    for (int iLine = 0; iLine < lineToFaceList.size(); iLine++)
    {
        LineToFace& lineToFace = lineToFaceList[iLine];
        //for boundary line, has only one face, the other is set to -1
        lineToFace.faceList[0] = -1;
        lineToFace.faceList[1] = -1;
    }
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        const Face& face = faceList[iFace];
        FaceToLine& faceToLine = faceToLineList[iFace];
        for (int iLine = 0; iLine < 3; iLine++)
        {
            Line lineLocal;
            for (int iPoint = 0; iPoint < 2; iPoint++)
            {
                //点按顺时针存取，lineLocal按逆序存储
                lineLocal.pointList[iPoint] = face.pointList[(iLine - iPoint + 2) % 3]; //ni xu
            }
            //log_info << "begin: findLine " << lineLocal.pointList[0] << " " << lineLocal.pointList[1];
            auto iter = findLine(lineMap, lineLocal);
            //log_info << "end: findLine";
            if (iter != lineMap.end())
            {
                //has found
                //log_info << "has found111";
                int iLineGlobal = iter->second;
                //log_info << "has found111 " << iLineGlobal;
                faceToLine.lineList[iLine] = iLineGlobal;
                //log_info << "has found222";
                if (lineToFaceList[iLineGlobal].faceList[0] == -1)
                {
                    lineToFaceList[iLineGlobal].faceList[0] = iFace;
                }
                else
                {
                    lineToFaceList[iLineGlobal].faceList[1] = iFace;
                }
                //log_info << "has found333";
                //log_info << "iFace = " << iFace << " iLineGlobal = " << iLineGlobal << " lineToFace[0 1] " << lineToFaceList[iLineGlobal].faceList[0] << " " << lineToFaceList[iLineGlobal].faceList[1];

            }
            else
            {
                //not found
                log_error << "lineToFace cal error";
            }
        }
    }


    //calculate faceToLineToFace from faceToLineList and lineToFaceList
    if (config->debugLevel > 1) log_info << "begin: faceToLineToFaceList";
    faceToLineToFaceList.resize(faceList.size());
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        FaceToLine& faceToLine = faceToLineList[iFace];
        FaceToLineToFace& faceToLineToFace = faceToLineToFaceList[iFace];
        for (int iLine = 0; iLine < 3; iLine++)
        {
            faceToLineToFace.lineList[iLine] = faceToLine.lineList[iLine];
            //log_info << "faceToLineToFace.lineList[iLine]: " << faceToLineToFace.lineList[iLine];
        }

        for (int iLine = 0; iLine < 3; iLine++)
        {
            int iLineGlobal = faceToLineToFace.lineList[iLine];
            LineToFace& lineToFace = lineToFaceList[iLineGlobal];
            if (lineToFace.faceList[0] == iFace)
            {
                if (lineToFace.faceList[1] == -1)
                {
                    //this face is boundary
                }
                else
                {
                    faceToLineToFace.faceList[iLine] = lineToFace.faceList[1];
                }

            }
            else
            {
                faceToLineToFace.faceList[iLine] = lineToFace.faceList[0];
            }

        }

    }

    //cal FaceToLineToFace other info: lineCenterPoint, normalList
    if (config->debugLevel > 1) log_info << "begin: cal FaceToLineToFace other info";
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        Face& face = faceList[iFace];
        Point& point0 = pointList[face.pointList[0]];
        Point& point1 = pointList[face.pointList[1]];
        Point& point2 = pointList[face.pointList[2]];


        //cal faceToLineToFace::centerPoint
        FaceToLineToFace& faceToLineToFace = faceToLineToFaceList[iFace];
        for (int i = 0; i < 3; i++)
        {
            faceToLineToFace.centerPoint.position(i) = 0.0;
        }
        for (int i = 0; i < 3; i++)
        {
            faceToLineToFace.centerPoint.position(i) += point0.position(i) / 3.0;
            faceToLineToFace.centerPoint.position(i) += point1.position(i) / 3.0;
            faceToLineToFace.centerPoint.position(i) += point2.position(i) / 3.0;
        }
        //log_info << "face center point position: " << iFace << " " << faceToLineToFace.centerPoint.position(0) << " " << faceToLineToFace.centerPoint.position(1) << " " << faceToLineToFace.centerPoint.position(2);

        //cal faceToLineToFace::lineCenterPointList
        for (int iLine = 0; iLine < 3; iLine++)
        {
            Point& lineCenterPoint = faceToLineToFace.lineCenterPointList[iLine];
            Line& line = lineList[faceToLineToFace.lineList[iLine]];
            for (int i = 0; i < 3; i++)
            {
                lineCenterPoint.position(i) = 0.0;
            }
            for (int i = 0; i < 3; i++)
            {
                for (int iPoint = 0; iPoint < 2; iPoint++)
                {
                    //log_info << line.pointList[iPoint];
                    Point& point = pointList[line.pointList[iPoint]];
                    lineCenterPoint.position(i) += point.position(i) * 0.5;
                    //log_info << " position " << point.position(i);
                }
                //log_info << "lineCenterPoint: " << lineCenterPoint.position(0) << " " << lineCenterPoint.position(1);
            }
            //log_info << "faceToLineToFace line center point position: " << faceToLineToFace.lineList[iLine] << " " << lineCenterPoint.position(0) << " " << lineCenterPoint.position(1) << " " << lineCenterPoint.position(2);
        }
        //exit(0);

        //cal faceToLineToFace::normalList(this line out normal)
        /*for (int iLine = 0; iLine < 3; iLine++)
        {
            Vector3d& normal = faceToLineToFace.normalList[iLine];
            Line& line = lineList[faceToLineToFace.lineList[iLine]];

            Point& point0 = pointList[line.pointList[0]];
            Point& point1 = pointList[line.pointList[1]];

            Eigen::Vector3d abNormal;
            Eigen::Vector3d faceCenterTolineCenter;

            faceCenterTolineCenter(0) = faceToLineToFace.lineCenterPointList[iLine].position(0) - faceToLineToFace.centerPoint.position(0);
            faceCenterTolineCenter(1) = faceToLineToFace.lineCenterPointList[iLine].position(1) - faceToLineToFace.centerPoint.position(1);
            faceCenterTolineCenter(2) = 0.0;

            abNormal(0) = point1.position(1) - point0.position(1);
            abNormal(1) = point0.position(0) - point1.position(0);
            abNormal(2) = 0.0;
            
            abNormal.normalize();

            if (faceCenterTolineCenter.dot(abNormal) > 0.0)
            {
                normal(0) = abNormal[0];
                normal(1) = abNormal[1];
                normal(2) = abNormal[2];
            }
            else
            {
                normal(0) = -abNormal[0];
                normal(1) = -abNormal[1];
                normal(2) = -abNormal[2];
            }
        }*/

    }
    //log_info << "end: cal FaceToLineToFace other info";


    //set lineCenterPoint for line in lineList
    if (config->debugLevel > 1) log_info << "begin: set lineCenterPoint for line in lineList";
    for (int iFace = 0; iFace < faceToLineToFaceList.size(); iFace++)
    {
        auto& faceToLineToFace = faceToLineToFaceList[iFace];
        for (int iLine = 0; iLine < 3; iLine++)
        {
            int iLineGlobal = faceToLineToFace.lineList[iLine];
            lineList[iLineGlobal].lineCenterPoint = faceToLineToFace.lineCenterPointList[iLine];

            //log_info << std::setprecision(14) << "lineList[iLineGlobal].lineCenterPoint: " << lineList[iLineGlobal].lineCenterPoint.position(0) << " " << lineList[iLineGlobal].lineCenterPoint.position(1);
        }
    }
    //======改到这里

    //==================================== begin: calculate for lineGroup =================================
    //确定lineGroup中单元在lineList中的序号
    //gmsh的只存了边界的line，lineGroup里边界的line序号也是边界line的，
    //EasyPIC里存的是所有line（包括边界和内部），下面是将lineGroup里line序号转换到全局line序号
    //line中点的顺序有六种可能，因此需要六次判断
    if (config->debugLevel > 1) log_info << "begin: set line index of lineGroup";
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        std::string lineGroupName = iter->first;
        auto& lineGroup = iter->second;
        lineGroup.lineList.resize(lineGroup.lineElementList.size());
        for (int i = 0; i < lineGroup.lineList.size(); i++)
        {
            int iLineLocal;
            Line line;
            line.pointList[0] = lineGroup.lineElementList[i]->pointList[0];
            line.pointList[1] = lineGroup.lineElementList[i]->pointList[1];

            auto iter = findLine(lineMap, line);
            if (iter != lineMap.end())
            {
                //has found
                int iLine = iter->second;
                lineGroup.lineList[i] = iLine;

                //log_info << "global lineGroup------------" << lineGroupName << " " << i << " " << iLine;
            }
            else
            {
                log_fatal << "not found, lineGroup: " << lineGroupName;
            }

        }
    }





    //calculate lineGroup of faceGroup's surface
    //if (config->debugLevel > 1) log_info << "begin: calculate lineGroup of faceGroup's surface";
    // !!!this process should be executed in geometry module of GUI
    //for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    //{
    //    string lineGroupName = iter->first + "Surface";
    //    auto& faceGroup = iter->second;

    //    vector<int> lineListTemp;

    //    map<int, int> faceMapTemp;
    //    for (int i = 0; i < faceGroup.faceList.size(); i++)
    //    {
    //        int iFace = faceGroup.faceList[i];
    //        faceMapTemp[iFace] = iFace;
    //    }

    //    for (int i = 0; i < faceGroup.faceList.size(); i++)
    //    {
    //        int iFace = faceGroup.faceList[i];
    //        FaceToLineToFace& faceToLineToFace = faceToLineToFaceList[iFace];

    //        for (int j = 0; j < 3; j++)
    //        {
    //            int iLine = faceToLineToFace.lineList[j];
    //            LineToFace& lineToFace = lineToFaceList[iLine];
    //            if (lineToFace.faceList[1] == -1)
    //            {
    //                continue;
    //            }

    //            int numIn = 0;
    //            int numWhich = -1;
    //            auto iter = faceMapTemp.find(lineToFace.faceList[0]);
    //            if (iter != faceMapTemp.end())
    //            {
    //                numIn++;
    //                numWhich = 0;
    //            }
    //            iter = faceMapTemp.find(lineToFace.faceList[1]);
    //            if (iter != faceMapTemp.end())
    //            {
    //                numIn++;
    //                numWhich = 1;
    //            }

    //            if (numIn == 1)
    //            {
    //                lineListTemp.push_back(iLine);
    //                lineList[iLine].normal = faceToLineToFace.normalList[j];
    //                lineList[iLine].isInGroup = true;

    //                //lineToFace.faceList[]的第一个face属于vacuum region
    //                if (numWhich == 0)
    //                {
    //                    int iFaceTemp = lineToFace.faceList[0];
    //                    lineToFace.faceList[0] = lineToFace.faceList[1];
    //                    lineToFace.faceList[1] = iFaceTemp;
    //                }
    //            }
    //        }
    //    }

    //    LineGroup lineGroup;
    //    lineGroup.name = lineGroupName;
    //    lineGroup.lineList = lineListTemp;
    //    lineGroupMap[lineGroupName] = lineGroup;
    //}


    //reorder faceToVolume.volumeList for objectSurface (the faceGroup corresponding to volumeGroup)
    //make sure faceToVolume.volumeList[0] is in vacuum region
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto& faceGroup = iter->second;

        std::map<int, int> faceMapTemp;
        for (int i = 0; i < faceGroup.faceList.size(); i++)
        {
            int iFace = faceGroup.faceList[i];
            faceMapTemp[iFace] = iFace;
        }

        for (int i = 0; i < faceGroup.faceList.size(); i++)
        {
            int iFace = faceGroup.faceList[i];
            FaceToLineToFace& faceToLineToFace = faceToLineToFaceList[iFace];

            for (int j = 0; j < 3; j++)
            {
                int iLine = faceToLineToFace.lineList[j];
                LineToFace& lineToFace = lineToFaceList[iLine];
                if (lineToFace.faceList[1] == -1)
                {
                    continue;
                }

                int numIn = 0;
                int numWhich = -1;
                auto iter = faceMapTemp.find(lineToFace.faceList[0]);
                if (iter != faceMapTemp.end())
                {
                    numIn++;
                    numWhich = 0;
                }
                iter = faceMapTemp.find(lineToFace.faceList[1]);
                if (iter != faceMapTemp.end())
                {
                    numIn++;
                    numWhich = 1;
                }

                if (numIn == 1)
                {
                    lineList[iLine].normal = faceToLineToFace.normalList[j];
                    lineList[iLine].isInGroup = true;

                    //lineToFace.faceList[]的第一个face属于vacuum region
                    if (numWhich == 0)
                    {
                        int iFaceTemp = lineToFace.faceList[0];
                        lineToFace.faceList[0] = lineToFace.faceList[1];
                        lineToFace.faceList[1] = iFaceTemp;
                    }
                }
            }
        }
    }



    //calculate pointList of lineGroup
    if (config->debugLevel > 1) log_info << "begin: calculate pointList of lineGroup";
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto& lineGroupName = iter->first;
        auto& lineGroup = iter->second;

        std::set<int> pointSet;
        for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
        {
            auto& line = lineList[lineGroup.lineList[iLine]];
            for (int iPoint = 0; iPoint < 2; iPoint++)
            {
                pointSet.insert(line.pointList[iPoint]);
            }
        }

        lineGroup.pointList.clear();
        for (auto iterPoint = pointSet.begin(); iterPoint != pointSet.end(); iterPoint++)
        {
            lineGroup.pointList.push_back(*iterPoint);
        }
    }

    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto& lineGroupName = iter->first;
        auto& lineGroup = iter->second;
        lineGroup.isPeriodic = false;

        lineGroup.numberLine = lineGroup.lineList.size();
        lineGroup.numberLineSerial = lineGroup.numberLine;
    }




    //check point number for periodic boundary lineGroup
    //reorder periodic lineGroup pointList
    if (config->debugLevel > 1) log_info << "begin: reorder periodic lineGroup pointList";
    for (int iLineGroup = 0; iLineGroup < config->mesh_param.mesh_periodic_edge0.size(); iLineGroup++)
    {
        std::string lineName0 = config->mesh_param.mesh_periodic_edge0[iLineGroup];
        std::string lineName1 = config->mesh_param.mesh_periodic_edge1[iLineGroup];
        auto& lineGroup0 = lineGroupMap.find(lineName0)->second;
        auto& lineGroup1 = lineGroupMap.find(lineName1)->second;

        //log_info << "lineName0 and lineName1: " << lineName0 << " " << lineGroup0.pointList.size() << " " << lineName1 << " " << lineGroup1.pointList.size();
        if (lineGroup0.pointList.size() != lineGroup1.pointList.size())
        {
            log_error << "periodic boundary lineGroup do not match: " << lineName0 << " " << lineGroup0.pointList.size() << " " << lineName1 << " " << lineGroup1.pointList.size();
            exit(0);
        }

        std::vector<int> pointListNew1;
        pointListNew1.resize(lineGroup0.pointList.size());
        for (int i = 0; i < lineGroup0.pointList.size(); i++)
        {
            int iPoint0 = lineGroup0.pointList[i];
            Point& point0 = pointList[iPoint0];

            double minLength = maxGeometryLength;
            int iPoint;
            for (int j = 0; j < lineGroup1.pointList.size(); j++)
            {
                int iPoint1 = lineGroup1.pointList[j];
                Point& point1 = pointList[iPoint1];
                double lineLength = calLineLength(point0, point1);
                if (lineLength < minLength)
                {
                    minLength = lineLength;
                    iPoint = iPoint1;
                }
            }
            pointListNew1[i] = iPoint;
            //log_info << "iPoint: " << iPoint;

        }
        lineGroup1.pointList = pointListNew1;

        //log_info << "pointListNew1.size: " << pointListNew1.size();
        //log_info << "lineName0 and lineName1: " << lineName0 << " " << lineGroup0.pointList.size() << " " << lineName1 << " " << lineGroup1.pointList.size();

        //log_info << "2222";
        //calculate periodicTranslationVector
        auto& point0 = pointList[lineGroup0.pointList[0]];
        auto& point1 = pointList[lineGroup1.pointList[0]];
        //log_info << "position: " << point1.position(0);

        //log_info << "2222 aaa";
        lineGroup0.periodicTranslationVector(0) = point1.position(0) - point0.position(0);
        //log_info << "2222 bbb";
        lineGroup0.periodicTranslationVector(1) = point1.position(1) - point0.position(1);
        lineGroup0.periodicTranslationVector(2) = point1.position(2) - point0.position(2);
        lineGroup1.periodicTranslationVector(0) = -lineGroup0.periodicTranslationVector(0);
        lineGroup1.periodicTranslationVector(1) = -lineGroup0.periodicTranslationVector(1);
        lineGroup1.periodicTranslationVector(2) = -lineGroup0.periodicTranslationVector(2);
        //log_info << "3333";
        lineGroup0.isPeriodic = true;
        lineGroup1.isPeriodic = true;

    }



    //reorder periodic lineGroup lineList
    if (config->debugLevel > 1) log_info << "begin: reorder periodic lineGroup lineList";
    for (int iLineGroup = 0; iLineGroup < config->mesh_param.mesh_periodic_edge0.size(); iLineGroup++)
    {
        std::string lineName0 = config->mesh_param.mesh_periodic_edge0[iLineGroup];
        std::string lineName1 = config->mesh_param.mesh_periodic_edge1[iLineGroup];
        auto& lineGroup0 = lineGroupMap.find(lineName0)->second;
        auto& lineGroup1 = lineGroupMap.find(lineName1)->second;
        //log_info << "++++++++++++++++++++++++++++++++++++++++++++++++++++++";
        //log_info << lineName0 << " " << lineName1;


        std::vector<int> lineListNew1;
        lineListNew1.resize(lineGroup0.lineList.size());

        
        for (int i = 0; i < lineGroup0.lineList.size(); i++)
        {
            int iLine0 = lineGroup0.lineList[i];
            Point& point0 = lineList[iLine0].lineCenterPoint;
            //log_info << "point0: " << point0.position(0) << " " << point0.position(1) << " " << point0.position(2);

            int iLine;
            double minLength = maxGeometryLength;
            for (int j = 0; j < lineGroup1.lineList.size(); j++)
            {
                int iLine1 = lineGroup1.lineList[j];
                Point& point1 = lineList[iLine1].lineCenterPoint;
                double lineLength = calLineLength(point0, point1);
                //log_info << "point1: " << point1.position(0) << " " << point1.position(1) << " " << point1.position(2) << " lineLength: " << lineLength;
                if (lineLength < minLength)
                {
                    minLength = lineLength;
                    iLine = iLine1;
                }
            }
            lineListNew1[i] = iLine;
            //log_info << "iLine0, iLine = " << iLine0 << " " << iLine << " minLength: " << minLength;
            Point& point1 = lineList[331].lineCenterPoint;
            //log_info << "point1: " << point1.position(0) << " " << point1.position(1) << " " << point1.position(2) << std::endl;


            lineList[iLine0].periodicLineIndex = iLine;
            lineList[iLine].periodicLineIndex = iLine0;

        }
        lineGroup1.lineList = lineListNew1;
    }
    //log_info << "end: reorder periodic lineGroup lineList";

    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        //slineGroup: lineGroup of serial mesh; pLineGroup: lineGroup of parallel mesh
        std::string sLineGroupName = iter->first;
        auto& sLineGroup = iter->second;

        for (int i = 0; i < sLineGroup.lineList.size(); i++)
        {
            int sLineNum = sLineGroup.lineList[i];
            //log_info << "Test serial mesh, ======, lineGroup-----  " << sLineGroupName << " " << sLineNum << " " << i;
        }
    }
    


    for (int iLine = 0; iLine < lineList.size(); iLine++)
    {
        lineList[iLine].isInGroup = false;
    }

    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto lineGroupName = iter->first;
        auto& lineGroup = iter->second;
        //log_info << lineGroupName << " line number: " << lineGroup.lineList.size();
        for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
        {
            //log_info << " ----------- line id: " << lineGroup.lineList[iLine];
            lineList[lineGroup.lineList[iLine]].isInGroup = true;
            lineList[lineGroup.lineList[iLine]].lineGroupName = lineGroupName;
        }

    }


    numBoundaryLine = 0;
    for (int iLine = 0; iLine < lineToFaceList.size(); iLine++)
    {
        auto& lineToFace = lineToFaceList[iLine];
        if (lineToFace.faceList[1] == -1)
        {
            numBoundaryLine++;
        }
    }
    if (config->debugLevel > 1) log_info << "end: calculate for lineGroup" << std::endl;
    //==================================== end: calculate for lineGroup =================================






    //cal area of line, 2d xy, the length of z direction is 1.0 m
    //2d rz, the length of phi direction is 2*pi*r m, r is r coordinate of line center point
    //if (config->debugLevel > 1) log_info << "cal area of line";
    //for (int iLine = 0; iLine < lineList.size(); iLine++)
    //{
    //    Line& line = lineList[iLine];
    //    Point& point0 = pointList[line.pointList[0]];
    //    Point& point1 = pointList[line.pointList[1]];

    //    Eigen::Vector3d vector0;
    //    vector0(0) = point1.position(0) - point0.position(0);
    //    vector0(1) = point1.position(1) - point0.position(1);
    //    vector0(2) = point1.position(2) - point0.position(2);

    //    if (config->isCylindricalCoordinateSystem)
    //    {
    //        line.area = vector0.norm() * 2.0 * config->const_pi * (point0.position(0) + point1.position(0)) / 2.0;
    //    }
    //    else
    //    {
    //        line.area = vector0.norm() * config->zDirectionLength;
    //    }
    //    
    //    //log_info << "line area: " << iLine << " " << line.area<< " " << line.pointList[0] << " " << line.pointList[1] << " " << line.pointList[2];
    //}


    ////cal volume of Face,  2d xy, the length of z direction is 1.0 m
    ////2d rz, the length of phi direction is 2*pi*r m, r is r coordinate of line center point
    //if (config->debugLevel > 1) log_info << "start calculate volume of face";
    //minVol = 1.0e99;
    //maxVol = 0.0;
    //totalVol = 0.0;
    //for (int iFace = 0; iFace < faceList.size(); iFace++)
    //{
    //    Face& face = faceList[iFace];
    //    Point& point0 = pointList[face.pointList[0]];
    //    Point& point1 = pointList[face.pointList[1]];
    //    Point& point2 = pointList[face.pointList[2]];

    //    face.vol = calTriangleArea2d(point0, point1, point2);

    //    totalVol += face.vol;

    //    if (config->isCylindricalCoordinateSystem)
    //    {
    //        face.vol *= 2.0 * config->const_pi * (point0.position(0) + point1.position(0) + point2.position(0)) / 3.0;
    //    }
    //    else
    //    {
    //        face.vol *= config->zDirectionLength;
    //    }

    //    if (face.vol > maxVol)
    //    {
    //        maxVol = face.vol;
    //    }

    //    if (minVol == 0.0)
    //    {
    //        minVol = face.vol;
    //    }
    //    if (face.vol < minVol)
    //    {
    //        minVol = face.vol;
    //    }

    //}
    averageVol = totalVol / faceList.size();





    //cal dualVolumeVolList
    //对于二维这里仍然是VolumeVol，而不是FaceArea，计算密度时候本质还是除以体积，我们规定二维直角坐标系，
    // 另一个方向的长度为 1.0 m
    //对于二维柱坐标系，另一个方向长度是2*pi*r
    if (config->debugLevel > 1) log_info << "start calculate dual volume";
    dualVolumeVolList.resize(pointList.size());
    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        dualVolumeVolList[iPoint] = 0.0;
    }
    for (int iFace = 0; iFace < faceList.size(); iFace++)
    {
        Face& face = faceList[iFace];
        Point& pointFaceCenter = faceToLineToFaceList[iFace].centerPoint;

        for (int iPoint = 0; iPoint < 3; iPoint++)
        {
            int iPointGlobal = face.pointList[iPoint];

            //point0: point for iPoint
            Point& point0 = pointList[iPointGlobal];
            Point& point1 = pointList[face.pointList[(iPoint + 1) % 3]];
            Point& point2 = pointList[face.pointList[(iPoint + 2) % 3]];

            Point pointLineCenter;

            //cal half trangle area for point1 side
            for (int i = 0; i < 3; i++)
            {
                pointLineCenter.position(i) = (point0.position(i) + point1.position(i)) / 2.0;
                pointLineCenter.position(i) = (point0.position(i) + point1.position(i)) / 2.0;
                pointLineCenter.position(i) = (point0.position(i) + point1.position(i)) / 2.0;
            }

            double volumeTemp = calTriangleArea2d(pointFaceCenter, point0, pointLineCenter);
            if (config->isCylindricalCoordinateSystem)
            {
                dualVolumeVolList[iPointGlobal] += volumeTemp * 2.0 * config->const_pi * (pointFaceCenter.position(0) + point0.position(0) + pointLineCenter.position(0)) / 3.0;
            }
            else
            {
                dualVolumeVolList[iPointGlobal] += volumeTemp * config->zDirectionLength;
            }
            

            //cal half trangle area for point2 side
            for (int i = 0; i < 3; i++)
            {
                pointLineCenter.position(i) = (point0.position(i) + point2.position(i)) / 2.0;
                pointLineCenter.position(i) = (point0.position(i) + point2.position(i)) / 2.0;
                pointLineCenter.position(i) = (point0.position(i) + point2.position(i)) / 2.0;
            }
            volumeTemp = calTriangleArea2d(pointFaceCenter, point0, pointLineCenter);
            if (config->isCylindricalCoordinateSystem)
            {
                dualVolumeVolList[iPointGlobal] += volumeTemp * 2.0 * config->const_pi * (pointFaceCenter.position(0) + point0.position(0) + pointLineCenter.position(0)) / 3.0;
            }
            else
            {
                dualVolumeVolList[iPointGlobal] += volumeTemp * config->zDirectionLength;
            }
            
        }
    }

    minDualVolumeVol = 1.0e99;
    maxDualVolumeVol = 0.0;

    for (int iPoint = 0; iPoint < pointList.size(); iPoint++)
    {
        double vol = dualVolumeVolList[iPoint];
        if (vol > maxDualVolumeVol)
        {
            maxDualVolumeVol = vol;
        }
        if (vol < minDualVolumeVol)
        {
            minDualVolumeVol = vol;
        }
    }



    ////calculate barycentric coordinate matrix 
    //// !!!后续删除串行网格的下面代码，省内存，因为粒子推动只用到了并行网格部分
    ////log_info << "begin calculate barycentric coordinate matrix ";
    //barycentric_coordinate_A_inv2d.resize(faceList.size());
    //TTensorList2d.resize(faceList.size());
    //detAList2d.resize(faceList.size());
    //for (int iFace = 0; iFace < faceList.size(); iFace++)
    //{
    //    Face& face = faceList[iFace];

    //    //eigen bug: inverse of MatrixXd is OK, but inverse of Matrix4d is wrong, but data store type of barycentric_coordinate_A_inv use Matrix4d, which 
    //    //should be OK.
    //    Eigen::MatrixXd matrix(3, 3);
    //    //MatrixXd matrix1(4,4);
    //    //Matrix4d matrix;
    //    for (int iPoint = 0; iPoint < 3; iPoint++)
    //    {
    //        Point& point = pointList[face.pointList[iPoint]];

    //        matrix(0, iPoint) = point.position(0);
    //        matrix(1, iPoint) = point.position(1);
    //        matrix(2, iPoint) = 1.0;
    //    }
    //    //matrix1 = matrix.inverse();
    //    barycentric_coordinate_A_inv2d[iFace] = matrix.inverse();


    //    Eigen::Matrix<double, 2, 2> matrixAbar;
    //    Eigen::Matrix<double, 2, 2> matrixTbar;
    //    Eigen::Matrix<double, 3, 2> matrixT;
    //    Eigen::Vector2d ab, ac;

    //    Eigen::Map<Eigen::Vector2d> a(pointList[face.pointList[0]].position);
    //    Eigen::Map<Eigen::Vector2d> b(pointList[face.pointList[1]].position);
    //    Eigen::Map<Eigen::Vector2d> c(pointList[face.pointList[2]].position);

    //    ab = b - a;
    //    ac = c - a;

    //    matrixAbar(0, 0) = ab(0);
    //    matrixAbar(1, 0) = ab(1);
    //    matrixAbar(0, 1) = ac(0);
    //    matrixAbar(1, 1) = ac(1);

    //    matrixTbar = abs(matrixAbar.determinant()) * matrixAbar.inverse();
    //    //matrixTbar = matrixAbar.determinant() * matrixAbar.inverse();

    //    matrixT(1, 0) = matrixTbar(0, 0);
    //    matrixT(1, 1) = matrixTbar(0, 1);
    //    matrixT(2, 0) = matrixTbar(1, 0);
    //    matrixT(2, 1) = matrixTbar(1, 1);


    //    matrixT(0, 0) = -(matrixTbar(0, 0) + matrixTbar(1, 0));
    //    matrixT(0, 1) = -(matrixTbar(0, 1) + matrixTbar(1, 1));

    //    TTensorList2d[iFace] = matrixT;

    //    detAList2d[iFace] = abs(matrixAbar.determinant());
    //    //detAList2d[iFace] = matrixAbar.determinant();

    //}
    //if (config->debugLevel > 1) log_info << "end calculate barycentric coordinate matrix ";




    //计算lineGroup或者faceGroup的总面积
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto& lineGroupName = iter->first;
        auto& lineGroup = iter->second;

        double totalAreaTemp = 0.0;
        for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
        {
            auto& line = lineList[lineGroup.lineList[iLine]];
            totalAreaTemp += line.area;
        }

        lineGroup.totalArea = totalAreaTemp;
    }



    outputMeshInfo();

    
}


void Mesh::outputMeshInfo()
{

    log_info << "========== serial Mesh Info ==============";
    log_info << " point number: " << pointList.size();
    log_info << " line number: " << lineList.size();
    log_info << " face number: " << faceList.size();
    log_info << " volume number: " << volumeList.size();
    log_info << "-------------------------------------";
    log_info << " dimension: " << config->dimension;
    log_info << " cell number: " << numberCell;
    log_info << "-------------------------------------";

    if (config->dimension == 3)
    {
        log_info << " boundary face number: " << numBoundaryFace;
    }
    else if (config->dimension == 2)
    {
        log_info << " boundary line number: " << numBoundaryLine;
    }
    
    
    log_info << " minVol: " << minVol;
    log_info << " maxVol: " << maxVol;
    log_info << " maxGeometryLength: " << maxGeometryLength;
    log_info << " minLineLength: " << minLineLength;
    log_info << " maxLineLength: " << maxLineLength;
    log_info << " minDualVolumeVol: " << minDualVolumeVol;
    log_info << " maxDualVolumeVol: " << maxDualVolumeVol;
    log_info << " geometry radius: " << radius;
    log_info << "---------------lineGroup----------------------";
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto lineGroupName = iter->first;
        auto& lineGroup = iter->second;
        log_info << lineGroupName << " lineGroup: " <<lineGroup.lineList.size();
    }
    log_info << std::endl;

    log_info << "---------------faceGroup----------------------";
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto faceGroupName = iter->first;
        auto& faceGroup = iter->second;
        log_info << faceGroupName << " faceGroup: " << faceGroup.faceList.size();
    }
    log_info << std::endl;

    log_info << "---------------volumeGroup----------------------";
    for (auto iter = volumeGroupMap.begin(); iter != volumeGroupMap.end(); iter++)
    {
        auto volumeGroupName = iter->first;
        auto& volumeGroup = iter->second;
        log_info << volumeGroupName << " volumeGroup: " << volumeGroup.volumeList.size();
    }
    log_info << std::endl;

    if (config->dimension == 3 && lineList.size() != 0)
    {
        log_warning << "3d line number should be 0";
    }
    else if (config->dimension == 2 && volumeList.size() != 0)
    {
        log_warning << "2d volume number should be 0";
    }

    if (minLineLength == 0.0 || minVol == 0.0 || minDualVolumeVol == 0.0)
    {
        log_error << "the mesh may have double point, double line, double face, or double volume, please check the mesh";
        exit(0);
    }
}


std::map<Line, int>::iterator Mesh::findLine(std::map<Line, int>& lineMap, Line& lineInp)
{
    std::map<Line, int>::iterator iterLine;
    Line line;

    line.pointList[0] = lineInp.pointList[0];
    line.pointList[1] = lineInp.pointList[1];
    iterLine = lineMap.find(line);
    if (iterLine != lineMap.end())
    {
        return iterLine;
    }

    line.pointList[0] = lineInp.pointList[1];
    line.pointList[1] = lineInp.pointList[0];
    iterLine = lineMap.find(line);
    if (iterLine != lineMap.end())
    {
        return iterLine;
    }

    return iterLine;
}



std::map<Face, int>::iterator Mesh::findFace(std::map<Face, int>& faceMap, Face& faceInp)
{
    std::map<Face, int>::iterator iterFace;
    Face face;

    face.pointList[0] = faceInp.pointList[0];
    face.pointList[1] = faceInp.pointList[1];
    face.pointList[2] = faceInp.pointList[2];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    face.pointList[0] = faceInp.pointList[0];
    face.pointList[1] = faceInp.pointList[2];
    face.pointList[2] = faceInp.pointList[1];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    face.pointList[0] = faceInp.pointList[1];
    face.pointList[1] = faceInp.pointList[0];
    face.pointList[2] = faceInp.pointList[2];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    face.pointList[0] = faceInp.pointList[1];
    face.pointList[1] = faceInp.pointList[2];
    face.pointList[2] = faceInp.pointList[0];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    face.pointList[0] = faceInp.pointList[2];
    face.pointList[1] = faceInp.pointList[0];
    face.pointList[2] = faceInp.pointList[1];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    face.pointList[0] = faceInp.pointList[2];
    face.pointList[1] = faceInp.pointList[1];
    face.pointList[2] = faceInp.pointList[0];
    iterFace = faceMap.find(face);
    if(iterFace != faceMap.end())
    {
        return iterFace;
    }

    return iterFace;
}



std::map<Volume, int>::iterator Mesh::findVolume(std::map<Volume, int>& volumeMap, Volume& volumeInp)
{
    std::map<Volume, int>::iterator iterVolume;
    Volume volume;

    std::vector<int> srcList;
    std::vector<std::vector<int>> desListList;
    srcList.resize(4);
    for (int i = 0; i < 4; i++)
    {
        srcList[i] = i;
    }
    perm(srcList, 0, 3, desListList);

    //log_info << "volume point list:  ";
    for (auto& desList : desListList)
    {
        //log_info << "volume point list:  " << desList[0] << " " << desList[1] << " " << desList[2] << " " << desList[3];

        volume.pointList[0] = volumeInp.pointList[desList[0]];
        volume.pointList[1] = volumeInp.pointList[desList[1]];
        volume.pointList[2] = volumeInp.pointList[desList[2]];
        volume.pointList[3] = volumeInp.pointList[desList[3]];
        iterVolume = volumeMap.find(volume);
        if (iterVolume != volumeMap.end())
        {
            return iterVolume;
        }
    }

    return iterVolume;
}

void Mesh::calBarycentricCoordinate(int iVolume, const Point& pointCurrent, std::vector<double>& barycentricCoordinate)
{
    /*barycentricCoordinate.resize(4);

    Eigen::Matrix4d& bc_coord_inv = barycentric_coordinate_A_inv[iVolume];
    Eigen::Vector4d x;
    Eigen::Vector4d bc;

    for (int i = 0; i < 3; i++)
    {
        x(i) = pointCurrent.position(i);
    }
    x(3) = 1.0;

    bc = bc_coord_inv * x;

    for (int i = 0; i < 4; i++)
    {
        barycentricCoordinate[i] = bc(i);
    }*/
}


void Mesh::calBarycentricCoordinate2d(int iVolume, const Point& pointCurrent, std::vector<double>& barycentricCoordinate)
{
    /*barycentricCoordinate.resize(3);

    Eigen::Matrix3d& bc_coord_inv = barycentric_coordinate_A_inv2d[iVolume];
    Eigen::Vector3d x;
    Eigen::Vector3d bc;

    for (int i = 0; i < 2; i++)
    {
        x(i) = pointCurrent.position(i);
    }
    x(2) = 1.0;

    bc = bc_coord_inv * x;

    for (int i = 0; i < 3; i++)
    {
        barycentricCoordinate[i] = bc(i);
    }*/
}


int Mesh::getVolumeIndex(const Point& point)
{
    int iVolumeCurrent = -1;

    int method = 0;

    //if (method == 0)
    //{
    //    //暴力循环，对每个volume进行判断，判断方法是质心坐标是否有负值。

    //    if (config->dimension == 3)
    //    {
    //        std::vector<double> bc(4);

    //        for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
    //        {
    //            calBarycentricCoordinate(iVolume, point, bc);
    //            if (bc[0] >= 0.0 && bc[1] >= 0.0 && bc[2] >= 0.0 && bc[3] >= 0.0)
    //            {
    //                iVolumeCurrent = iVolume;

    //                auto& pointCenter = volumeToFaceToVolumeList[iVolume].centerPoint;
    //                log_info << "point position: " << pointCenter.position(0) << " " << pointCenter.position(1) << " " << pointCenter.position(2);

    //                break;
    //            }
    //        }
    //    }
    //    else if (config->dimension == 2)
    //    {
    //        std::vector<double> bc(3);

    //        for (int iFace = 0; iFace < faceList.size(); iFace++)
    //        {
    //            calBarycentricCoordinate2d(iFace, point, bc);
    //            if (bc[0] >= 0.0 && bc[1] >= 0.0 && bc[2] >= 0.0)
    //            {
    //                iVolumeCurrent = iFace;

    //                auto& pointCenter = faceToLineToFaceList[iFace].centerPoint;
    //                log_info << "point position: " << pointCenter.position(0) << " " << pointCenter.position(1);

    //                break;
    //            }
    //        }
    //    }
    //    else
    //    {
    //        log_error << "the dimension not support: " << config->dimension;
    //        exit(0);
    //    }


    //}
    //else if (method == 1)
    //{
    //    //下面代码只能追踪模拟区域内部没有几何实体的，因为是追踪的直线，如果有几何实体，直线可能达到几何实体表面，而找不到目标点坐在volume
    ////============================= method 1========================================================================================
    //    Point pointBegin = volumeToFaceToVolumeList[0].centerPoint;
    //    Point pointCurrent = pointBegin;

    //    std::vector<double> displacement(3);
    //    displacement[0] = point.position(0) - pointCurrent.position(0);
    //    displacement[1] = point.position(1) - pointCurrent.position(1);
    //    displacement[2] = point.position(2) - pointCurrent.position(2);

    //    std::vector<double> barycentricCoordinate(4);
    //    calBarycentricCoordinate(iVolumeCurrent, pointCurrent, barycentricCoordinate);

    //    double fraction = 0.0;




    //    // equation(12)
    //    Eigen::Vector3d Txd;
    //    Eigen::Vector4d Tx1;
    //    Eigen::Vector4d u;

    //    while (true)
    //    {
    //        if (fraction == 1.0)
    //        {
    //            break;
    //        }

    //        if (displacement[0] * displacement[0] + displacement[1] * displacement[1] + displacement[2] * displacement[2] == 0.0)
    //        {
    //            //log_info << "displacement[0] == 0.0";
    //            fraction = 1.0;
    //            break;
    //        }

    //        auto& TTensor = TTensorList[iVolumeCurrent];
    //        double detA = detAList[iVolumeCurrent];

    //        VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolumeCurrent];


    //        Eigen::Map<Eigen::Vector3d> xd(displacement.data());
    //        Eigen::Map<Eigen::Vector4d> y(barycentricCoordinate.data());
    //        Tx1 = TTensor * xd;
    //        //cout << "------------" <<TTensor <<xd;
    //        //cout << "------------" <<xd;
    //        //cout << "------------" <<u;

    //        double uMin = std::isnormal(detA) && detA <= 0.0 ? 1.0e99 : 1.0 / detA;
    //        int i_face_cross = -1;
    //        for (int i = 0; i < 4; i++)
    //        {
    //            if (Tx1(i) < -detA * 1.0e-15)
    //            {
    //                u(i) = -y(i) / Tx1(i);
    //                if (0.0 <= u(i) && u(i) < uMin)
    //                {
    //                    uMin = u(i);
    //                    i_face_cross = i;
    //                }
    //            }

    //        }


    //        //move particle in barycentric coordinate
    //        y = y + uMin * Tx1;

    //        //clamp to zero any negative coordiantes generated by round-off error
    //        for (int i = 0; i < 4; i++)
    //        {
    //            if (i == i_face_cross)
    //            {
    //                y(i) = 0.0;
    //            }
    //            else
    //            {
    //                if (y(i) <= 0.0)
    //                {
    //                    y(i) = 0.0;
    //                }
    //            }
    //        }

    //        if (i_face_cross == -1)
    //        {
    //            //log_info << "i_face_cross = -1";
    //            y = y / y.sum();

    //        }

    //        //log_info << "---------------------------" ;
    //        //log_info << "uMin: " << uMin << " detA = " << detA;
    //        if (i_face_cross == -1)
    //        {
    //            //particle stay in the cell
    //            fraction = 1.0;
    //        }
    //        else
    //        {
    //            double f = detA * uMin;
    //            displacement[0] *= (1.0 - f);
    //            displacement[1] *= (1.0 - f);
    //            displacement[2] *= (1.0 - f);

    //            fraction = fraction + f * (1.0 - fraction);

    //            FaceToVolume& faceToVolume = faceToVolumeList[volumeToFaceToVolume.faceList[i_face_cross]];
    //            int iFaceGlobal = volumeToFaceToVolume.faceList[i_face_cross];
    //            Face& face = faceList[iFaceGlobal];

    //            if (face.isInGroup)
    //            {
    //                break;
    //            }
    //            else
    //            {
    //                int iVolumeNeighbour = volumeToFaceToVolume.volumeList[i_face_cross];
    //                //log_info << "calCartesianCoordinateOfParticle";
    //                pointCurrent.position(0) = pointBegin.position(0) + fraction * (point.position(0) - pointBegin.position(0));
    //                pointCurrent.position(1) = pointBegin.position(1) + fraction * (point.position(1) - pointBegin.position(1));
    //                pointCurrent.position(2) = pointBegin.position(2) + fraction * (point.position(2) - pointBegin.position(2));

    //                iVolumeCurrent = iVolumeNeighbour;
    //                calBarycentricCoordinate(iVolumeCurrent, pointCurrent, barycentricCoordinate);

    //            }
    //        }

    //    }

    //    //============================= method 1 end ========================================================================================
    //}
    

    return iVolumeCurrent;
}

int Mesh::getLocalFaceIndex(int iVolume, int iFaceGlobal)
{
    VolumeToFaceToVolume& volumeToFaceToVolume = volumeToFaceToVolumeList[iVolume];
    for (int i = 0; i < 4; i++)
    {
        //log_info << "-----------------" << i;
        //log_info << "iFaceGlobal: " << iFaceGlobal << " volumeToFaceToVolume.faceList[i]: " << volumeToFaceToVolume.faceList[i];
        if (volumeToFaceToVolume.faceList[i] == iFaceGlobal)
        {
            return i;
        }
    }
    return -1;
}

int Mesh::getLocalLineIndex(int iFace, int iLineGlobal)
{
    FaceToLineToFace& faceToLineToFace = faceToLineToFaceList[iFace];
    for (int i = 0; i < 3; i++)
    {
        //log_info << "-----------------" << i;
        //log_info << "iFaceGlobal: " << iFaceGlobal << " volumeToFaceToVolume.faceList[i]: " << volumeToFaceToVolume.faceList[i];
        if (faceToLineToFace.lineList[i] == iLineGlobal)
        {
            return i;
        }
    }
    return -1;
}

int Mesh::getLocalPointIndexInVolume(const int& iVolume, const int& iPointGlobal)
{
    Volume& volume = volumeList[iVolume];
    for (int i = 0; i < 4; i++)
    {
        //log_info << "-----------------" << i;
        //log_info << "iFaceGlobal: " << iFaceGlobal << " volumeToFaceToVolume.faceList[i]: " << volumeToFaceToVolume.faceList[i];
        if (volume.pointList[i] == iPointGlobal)
        {
            return i;
        }
    }
    return -1;
}


int Mesh::getLocalPointIndexInFace(const int& iFace, const int& iPointGlobal)
{
    Face& face = faceList[iFace];
    for (int i = 0; i < 3; i++)
    {
        //log_info << "-----------------" << i;
        //log_info << "iFaceGlobal: " << iFaceGlobal << " volumeToFaceToVolume.faceList[i]: " << volumeToFaceToVolume.faceList[i];
        if (face.pointList[i] == iPointGlobal)
        {
            return i;
        }
    }
    return -1;
}


void Mesh::setInternalFaceGroup()
{
    for (auto iter = faceGroupMap.begin(); iter != faceGroupMap.end(); iter++)
    {
        auto faceGroupName = iter->first;
        auto& faceGroup = iter->second;

        faceGroup.isInternal = false;

        std::vector<int> neighborVolumeNumberList;
        std::vector<Face> faceListTemp;
        std::map<Face, int> faceMapTemp;

        auto& faceElementList = faceGroup.faceElementList;
        neighborVolumeNumberList.resize(faceElementList.size());

        for (int iFaceElement = 0; iFaceElement < faceElementList.size(); iFaceElement++)
        {
            neighborVolumeNumberList[iFaceElement] = 0;

            Element* faceElement = faceElementList[iFaceElement];
            Face face;
            face.pointList[0] = faceElement->pointList[0];
            face.pointList[1] = faceElement->pointList[1];
            face.pointList[2] = faceElement->pointList[2];
            faceListTemp.push_back(face);

            faceMapTemp[face] = iFaceElement;
        }
        

        for (int iVolume = 0; iVolume < volumeList.size(); iVolume++)
        {
            const Volume& volume = volumeList[iVolume];
            for (int iFace = 0; iFace < 4; iFace++)
            {
                Face faceLocal;
                for (int iPoint = 0; iPoint < 3; iPoint++)
                {
                    //点按顺时针存取，faceLocal按逆序存储
                    //faceLocal.pointList[iPoint] = volume.pointList[(iFace+iPoint+1)%4] //shun xu
                    faceLocal.pointList[iPoint] = volume.pointList[(iFace - iPoint + 3) % 4]; //ni xu
                }
                auto iter = findFace(faceMapTemp, faceLocal);
                if (iter != faceMapTemp.end())
                {
                    //has found
                    int iFaceGlobal = iter->second;
                    neighborVolumeNumberList[iFaceGlobal]++;

                }
                else
                {
                    //not found
                    
                }

            }

        }

        for (int i = 0; i < neighborVolumeNumberList.size(); i++)
        {
            //log_info << faceGroupName << " " << neighborVolumeNumberList[i];
            if (neighborVolumeNumberList[i] == 2)
            {
                faceGroup.isInternal = true;
                if (config->debugLevel > 1) log_info << "faceGroup " << faceGroupName << " is internal";
                break;
            }
        }


    }

    
}

void Mesh::setInternalLineGroup()
{
    for (auto iter = lineGroupMap.begin(); iter != lineGroupMap.end(); iter++)
    {
        auto lineGroupName = iter->first;
        auto& lineGroup = iter->second;

        lineGroup.isInternal = false;

        std::vector<int> neighborFaceNumberList;
        std::vector<Line> lineListTemp;
        std::map<Line, int> lineMapTemp;

        auto& lineElementList = lineGroup.lineElementList;
        neighborFaceNumberList.resize(lineElementList.size());

        for (int iLineElement = 0; iLineElement < lineElementList.size(); iLineElement++)
        {
            neighborFaceNumberList[iLineElement] = 0;

            Element* lineElement = lineElementList[iLineElement];
            Line line;
            line.pointList[0] = lineElement->pointList[0];
            line.pointList[1] = lineElement->pointList[1];
            lineListTemp.push_back(line);

            lineMapTemp[line] = iLineElement;
        }


        for (int iFace = 0; iFace < faceList.size(); iFace++)
        {
            const Face& face = faceList[iFace];
            for (int iLine = 0; iLine < 3; iLine++)
            {
                Line lineLocal;
                for (int iPoint = 0; iPoint < 2; iPoint++)
                {
                    //点按顺时针存取，lineLocal按逆序存储
                    //faceLocal.pointList[iPoint] = volume.pointList[(iFace+iPoint+1)%4] //shun xu
                    lineLocal.pointList[iPoint] = face.pointList[(iLine - iPoint + 2) % 3]; //ni xu
                }
                auto iter = findLine(lineMapTemp, lineLocal);
                if (iter != lineMapTemp.end())
                {
                    //has found
                    int iLineGlobal = iter->second;
                    neighborFaceNumberList[iLineGlobal]++;

                }
                else
                {
                    //not found

                }

            }

        }

        for (int i = 0; i < neighborFaceNumberList.size(); i++)
        {
            //log_info << faceGroupName << " " << neighborVolumeNumberList[i];
            if (neighborFaceNumberList[i] == 2)
            {
                lineGroup.isInternal = true;
                if (config->debugLevel > 1) log_info << "lineGroup " << lineGroupName << " is internal";
                break;
            }
        }


    }
}


bool Point::operator < (const Point& point) const
{
    double small = lengthSmallValue;
    if (position(0) + small < point.position(0))
    {
        return true;
    }
    else if (fabs(position(0) - point.position(0)) < small)    //(position(0) == point.position(0))
    {
        if (position(1) + small < point.position(1))    //position(1) < point.position(1)
        {
            return true;
        }
        else if (fabs(position(1) - point.position(1)) < small)            // position(1) == point.position(1)
        {
            if (position(2) + small < point.position(2))       //position(2) < point.position(2)
            {
                return true;
            }
            else
            {
                return false;
            }
        }
        else
        {
            return false;
        }
    }
    else
    {
        return false;
    }
}


VolumeInCylindricalCoordSysCalculator::VolumeInCylindricalCoordSysCalculator()
{
    IntRules = new mfem::IntegrationRules(0, mfem::Quadrature1D::GaussLegendre);
    ir = &IntRules->Get(mfem::Geometry::TRIANGLE, 10);
    Tr = new mfem::IsoparametricTransformation();
}

void VolumeInCylindricalCoordSysCalculator::setPoints(Point& p1, Point& p2, Point& p3)
{

    log_info << "aaa";
    //set Tr pointMat
    mfem::DenseMatrix& pointMat = Tr->GetPointMat();
    pointMat.SetSize(3, 3);
    log_info << "bbb";
    pointMat(0, 0) = p1.position(0);
    pointMat(1, 0) = p1.position(1);
    pointMat(2, 0) = 0.0;

    pointMat(0, 1) = p2.position(0);
    pointMat(1, 1) = p2.position(1);
    pointMat(2, 1) = 0.0;

    pointMat(0, 2) = p3.position(0);
    pointMat(1, 2) = p3.position(1);
    pointMat(2, 2) = 0.0;
    log_info << "ccc";
    //set Tr FE
    Tr->SetFE(&TriangleFE);
    log_info << "ddd";

}

//ref mfem void DomainLFIntegrator::AssembleRHSElementVect(const FiniteElement &el,ElementTransformation& Tr,Vector& elvect)
double VolumeInCylindricalCoordSysCalculator::cal()
{
    double area = 0.0;
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
        log_info << "ip: " << i;
        const mfem::IntegrationPoint& ip = ir->IntPoint(i);

        Tr->SetIntPoint(&ip);

        //get integration point coordinate
        double x[3];
        mfem::Vector transip(x, 3);
        Tr->Transform(ip, transip);
        log_info << "ip coordinate: " << x[0] << " " << x[1] << " " << x[2];

        double val = Tr->Weight() * (x[0]) * 2.0 * 3.1415926;

        area += val * ip.weight;
    }
    return area;
}




